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HYDROGEN  STARK  BROADENING  CALCULATIONS  WITH 
THE  UNIFIED  CLASSICAL  PATH  THEORY 
C.  R.  Vidal,  J.  Cooper,  and  E.  W,  Smith 


The  unified  theory  has  been  generalized  for  the 
case  of  upper  and  lower  state  interaction  by  introducing 
a  more  compact  tetradic  notation.    The  general  result 
is  then  applied  to  the  Stark  broadening  of  hydrogen. 
The  thermal  average  of  the  time  development  operator 
for  upper  and  lower  state  interaction  is  presented.  Ex- 
cept for  the  time  ordering  it  contains  the  effect  of  finite 
interaction  time  between  the  radiator  and  perturbers  to 
all  orders,  thus  avoiding  a  Lewis  type  cutoff.    A  simple 
technique  for  evaluating  the  Fourier  transform  of  the 
thermal  average  has  been  developed.    The  final  calcu- 
lations based  on  the  unified  theory  and  on  the  one-elect- 
ron theory  are  compared  with  measurements  in  the  high 
and  low  electron  density  regime.     The  unified  theory  cal- 
culations cover  the  entire  line  profile  from  the  line  center 
to  the  static  wing  and  the  simpler  one-electron  theory 
calculations  provide  the  line  intensities  only  in  the  line 
wings. 

Key  words:    Classical  path;  hydrogen  lines;  line 
wings;  one-electron  theory;  Stark 
broadening;  unified  theory. 
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I.  INTRODUCTION 

For  the  first  few  Balmer  lines  of  hydrogen,  recent  papers 

(Cerardoand  Hill,  1966;  Bacon  and  Edwards,   1968;  Kepple  and  Griem, 

1968;  Birkeland,  Oss  and  Braun,  1969)  have  demonstrated  fairly  good 

agreement  between  measurements  in  high  electron  density  plasmas 
16-3 

(n   >  10     cm     )  and  improved  calculations  of  the  so  called  "modified 
e 

impact  theory".     The  experimental  and  theoretical  half-widths  differ 

less  than  about  10%.    However,  measurements  of  the  Lyman-a  wings 

(Boldt  and  W.  Cooper,  1964;  Elton  and  Griem,  1964)  and  low  electron 

13  -3 

density  measurements  (n   =10     cm     )  of  the  higher  Balmer  and 

e 

Paschen  lines  (Ferguson  and  Schliiter,  1963;  Vidal,  1964;  Vidal,  1965) 
have  revealed  parts  of  the  hydrogen  line  profile,  for  which  the  modi- 
fied impact  theory  appears  to  break  down.    For  the  higher  series 
members  better  agreement  has  been  obtained  with  quasi- static  cal- 
culations (Vidal,   1965).     The  reason  the  current  impact  theories 
break  down  is  that  these  theories  correct  the  completed  collision 
assumption  by  means  of  the  Lewis  cutoff  (Lewis,   1961)  which  is  only 
correct  to  second  order.    With  this  cutoff  it  was  possible  to  extend 
the  range  of  validity  for  the  impact  theory  beyond  the  plasma  frequency. 
However,   in  the  distant  wings,  where  the  electron  broadening  becomes 
quasistatic,  the  second  order  perturbation  treatment  with  the  Lewis 
cutoff  breaks  down  because  the  time  development  operator  must  then 
be  evaluated  to  all  orders.    Attempts  to  correct  the  second  order  theory 
have  been  made  already  (Griem,   1965;  Shen  and  J.  Cooper,  1969),  but 
these  theories  still  make  the  completed  collision  assumption  by  re- 
placing the  time  development  operator  by  the  corresponding  S-matrix, 
and  so  it  has  to  be  emphasized  that  in  conjunction  with  the  Lewis  cutoff 
these  theories  would  only  be  correct  to  second  order.     The  impact 
theory  in  its  present  form  is  intrinsicly  not  able  to  describe  the  static 
wing  and  the  transition  region  to  the  line  center  where  dynamic  effects 
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'  cannot  be  neglected.     To  overcome  this  problem,   several  semiemp- 
irical  procedures  (Griem,  1962;  Griem,  1967a;  F.  Edmonds,  Schluter 
and  Wells,   1967)  have  been  suggested  to  generate  a  smooth  transition 
from  the  modified  impact  theory  to  the  static  wing. 


been  reinvestigated  in  two  review  papers  (  E.  Smith,   Vidal  and  J. 
Cooper,   1969a,   1969b),  which  are  from  now  on  referred  to  as  papers 
I  and  II.     The  purpose  of  I  and  II  was  to  state  clearly  the  different 
approximations  which  are  required  to  obtain  the  classical  path  theories 
of  line  broadening  and  to  find  out  where  these  theories  are  susceptible 
to  improvements.    In  a  manner  similar  to  the  Mozer-Baranger  treat- 
ment of  electric  microfield  distribution  functions  (Baranger  and  Mozer, 
1959,  I960),  it  was  shown  that  the  general  thermal  average  can  be 
expanded  in  two  ways,  one  of  which  leads  to  the  familiar  impact  tieory 
describing  the  line  center  (Baranger,  1958,  1962;  Griem,  Kolb  and 
Shen,   1959,   1962).     The  other  expansion  represents  a  generalized 
version  of  the  one  electron  theory  (J.   Cooper,   1966),  which  holds  in 
the  line  wings.    It  is  also  shown  that  there  is  generally  a  considerable 
domain  of  overlap  between  the  modified  impact  theory  and  the  one 
electron  theory.    Based  on  these  results,  a  "unified  theory"  was  then 
developed  (E.  Smith,  J.  Cooper  and  Vidal,   1969),  henceforth  referred 
to  as  paper  III,  which  presents  the  first  line  shape  expression  which  is 
valid  from  the  line  center  out  to  the  static  line  wing  including  the 
problematic  transition  region.     The  line  shape  obtained  by  the  unified 
theory  has  the  form 


Recently  the  classical  path  methods  in  line  broadening  have 


(I.  1) 
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where  d,  Auu  and  jb(Aaj)  are  operators.    In  paper  III  it  was  shown  that 
the  familiar  impact  theories,  which  hold  in  the  line  center,  may  be 
obtained  by  making  a  Markoff  approximation  in  the  unified  theory, 
while  the  one  electron  theory  describing  the  line  wings  is  just  a  wing 
expansion  of  the  unified  theory.     Consequently  the  crucial  problem  for 
any  line  broadening  calculation  is  to  evaluate  the  matrix  elements  of 
jL(AtD),  which  is  essentially  the  Fourier  transform  of  the  thermal  average 
(see  Eq.  (46)  and  (47)  of  paper  III).     This  will  be  done  in  detail  in  this 
paper  for  the  general  case  of  upper  and  lower  state  interactions. 

In  the  following  Sec.  II  we  start  with  a  brief  summary  of  the 

basic  relations  which  are  required  for  the  classical  path  approach 

pursued  here.    We  then  generalize  the  results  of  the  unified  theory  to 

include  lower  state  interaction  (Sec.  IV)  after  introducing  a  more 

compact  tetradic  notation  (Sec.  III).    From  this  general  result  we  turn 

to  the  specific  problem  of  hydrogen  by  discussing  briefly  the  no 

quenching  assumption  (Sec,  V)  and  deriving  the  thermal  average  5^^^(t) 

(see  Eq.  (47)  of  paper  III)  for  the  general  case  of  upper  and  lower  state 

interaction  (Sec.  VI).    We  next  investigate  the  multipole  expansion  of 

the  classical  interaction  potential  in  the  time  development  operator 

(1 ) 

(Sec.  VII).     The  thermal  average  5      (t)  is  then  evaluated  in  two  steps 

by  first  performing  a  spherical  average  (Sec.  VIII)  and  then  an  average 

over  the  collision  parameters:  some  reference  time  t  ,   impact  para- 

o 

meter  p  and  velocity  v  (Sec.  IX).     Appendix  A  gives  the  computer 
program  which  we  used  in  calculating  the  thermal  average  for  dipole 
interactions  including  lower  state  interactions.     The  large  time  limit 
of  the  thermal  average,  which  leads  to  the  familiar  impact  theories  in 
the  line  center,  is  investigated  in  detail  in  Appendix  B  for  different 
cutoff  procedures  and  compared  with  the  results  in  the  literature.  In 
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Sec.  X,  a  method  is  developed  for  performing  the  Fourier  transform- 
ation of  the  thermal  average  and  it  leads  us  to  the  crucial  function  for 
any  classical  path  theory  of  Stark  broadening.    This  function  is  finally 
applied  in  Sec.  XI  to  the  one  electron  theory,  which  forms  the  basis 
for  the  asymptotic  wing  expansion,  and  in  Sec.  XII  to  the  unified  theory, 
which  describes  the  whole  line  profile  from  the  line  center  to  the  static 
wing.    Numerical  results  are  given  for  the  hydrogen  line  profiles  as 
measured  by  Boldt  and  W.  Cooper,  1964;  Elton  and  Griem,  1964,  and 
j    Vidal,  1964,  1965.     The  computer  program  for  the  unified  theory  cal- 
culations and  the  asymptatic  wing  expansion  is  given  and  explained  in 
Appendix  C. 

11.    BASIC  RELATIONS 

In  this  section  we  will  briefly  outline  the  basic  relations  which 
^  are  used  in  our  classical  path  treatment  of  line  broadening. 

As  discussed  in  Sec,  2  of  paper  I,  we  are  considering  a  system 
containing  a  single  radiator  and  a  gas  of  electrons  and  ions.    We  will 
make  the  usual  quasi- static  approximation  for  the  ions  by  regarding 
their  electric  field  6.  as  being  constant  during  the  times  of  interest 
1/Au)  .    This  approximation  is  usually  very  good  because  the  region 
where  ion  dynamics  are  important  is  normally  well  inside  the  half 
width  of  the  line  except  for  a  few  cases  such  as  the  n-a  lines  of 
hydrogen  (Griem,  1967b).     The  complete  line  profile  I(uu )  is  then  given 
by  the  microfield  average  (see  Eq.  (3)  of  paper  II) 

00 

i((w)=  /'p(e.)  I  ((ju,e.)  de.  (II.  1) 

where  the  normalized  distribution  function  P{F..)  is  the  low  frequency 

component  of  the  fluctuating  electric  microfields.  Due  to  shielding  effects 

Pit.)  depends  on  the  shielding  parameter  r  /D  where  r    and  D   are  the 
1  o  o 

mean  particle  distance  and  the  Debye  length  (for  electrons  CTnl^jrespectively. 
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With  the  static  ion  approximation  we  have  reduced  the  problem 

to  a  calculation  of  the  electron  broadening  of  a  radiator  in  a  static 

electric  fields  .     The  resulting  line  profile  I(aj,e.)  is  then  simply 
1  1 

averaged  over  all  possible  ion  fields  to  give  the  complete  line  profile 

I((jd).     The  static  ion  field  will  be  used  to  define  the  z-axis  for  the 

radiator  and  the  ion- radiator  interaction  will  be  taken  to  be  the  dipole 

interaction  eZP, .  where  -eZ  denotes  the  Z-component  of  the  radiators 
1 

dipole  moment. 

If  the  unperturbed  radiator  is  described  by  a  Hamiltonian  H  , 

a 

we  may  then  define  a  Hamiltonian  for  a  radiator  in  the  static  field 


P..  by 


H    =  H    +  eZp,.  {II.  2) 

o        a  1 


The  complete  Hamiltonian  for  the  system  is  then  given  by 


H  =  H    +  V  (R,  X,  V,  t)  (II.  3) 

o  e 


where        denotes  the  electron  radiator  interaction.    In  this  equation, 
X  and  V  are  3N  vectors  x  =  (x^ ,  5c^,   •  .  .  ,  x^),  v  =  (v^,^^'   '  *  *  '  ^n^' 
which  denote  the  positions  and  velocities  of  the  N  electrons  and  R 
denotes  some  internal  radiator  coordinates.    For  one-electron  atoms, 
R  is  the  position  of  the  "orbital"  electron  relative  to  the  nucleus. 
The  interaction  V    will  be  regarded  as  a  sum  of  binary  interactions, 

V    (R,  X,  V,  t)=   y  V.(R,  X.,  v.,  t)  (II.  4) 

J 

where        denotes  the  interaction  between  the  radiator  and  a  single 

electron.    As  is  well  known  the  line  shape  I(aJ,£.)  may  be  given  by 

the  Fourier  transform  of  an  autocorrelation  function  C(t)  (Baranger,  1962) 

fe'^^C(t)dt  (II.  5) 

o 


In  the  classical  path  approximation,  the  correlation  function  for  electric 
dipole  radiation  is  given  by 

C(t)=Tr    I  d  (t"'"  (t)  d  T  (t)  >     p     \    ,  (II.  6) 

a  I         a  a       av    a  / 

where       and   p    denote  the  dipole  moment  and  the  density  matrix  for 
the  radiator.     The  thermal  average  denoted  by  the  subscript  av 
represents  the  average  over  electron  states  (see  Eq.  (47)  of  paper  I): 


<T"'"(t)dT  (t))      =/dxdvP  (x)W(v)  t'^(R,x,  V,  t)  d  T  (R,  x,  v,  t)      (II.  7) 
a  a       av    I  a  a 


where  P(x)  and  W(v)  are  the  position  and  velocity  distribution  functions 
for  the  electron  perturbers  (defined  by  Eqs.(37)  to    (40)  in  paper  II). 
The  time  development  operator  for  the  system  T  (R,  x,  v,  t)  is  the 
solution  of  the  differential  equation 

ih^T  (t)=   [h    +V  (t)l    T  (t)  (11.  8) 

9t    a  L     °        e      I  a 

and  it  may  be  written  in  an  interaction  representation  defined  by 


T  (R,  X,  V,  t)  =  exp  (-itH  /h)  U  (R,  x,  v,  t)  (II.  9) 

a  o.  a 


\  where 


fti-^  u  (t)  =  V  (t)  u  (t)  (11.10) 

3t      a  e  a 
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and 


V  (t)  =  exp  (itH  /h)V  (t)  exp  (-itH  /h)  .         (II.  11) 

e  o        e  o 


It  should  be  noted  that  V  (t)  is  identical  with  V  (t)  in  paper  II  except 

e  e 

that  we  have  not  yet  made  the  no  quenching  assumption  which  removes 

the  unperturbed  part  H    in  the  Hamiltonian  H    in  Eq.  (II.  11),  Using 

a  o 

the  time  ordering  operator  O,  U  (t)  may  be  written  in  the  form 

U    (R,  X, 'v,  t)  =      exp  f  -  ^  /'v  (R,  X,  V,  t')  dt'   !.  (11.12) 


a 


To  evaluate  the  trace  over  atomic  states  in  Eq.  (II.  6),   it  is  convenient 

to  use  the  H    eigenstates  I  a 
o 

Hence,  using  U  (t)  we  have 


to  use  the  H    eigenstates  la),    lb),   •  •  •  with  the  eigenvalues  E  ,  E  , 
o  °  a  b 


-iju  ,  t 
dc 


C(t)=  ^<a!d|b)  (crd|d>e  (11.13) 
abed 


(b|u"^(t)lc)  <dl  U  (t)|a> 

2i  Si 


<a|p  la) 
av  a 


where 


0)  ,    =  (E,-E  )/h  (II.  14) 

dc  d  c 


In  paper  II  and  III,  the  correlation  function  C(t)  was  evaluated 
for  the  case  of  no  lower  state  interactions  in  order  to  keep  the  math- 
ematics as  simple  as  possible  because  one  of  the  U  (t)  operators  in 

3. 

Eq.  (II.  13)  may  then  be  replaced  by  a  unit  operator.    In  this  paper  we 
will  give  a  more  general  evaluation  of  C(t)  which  includes  lower  state 
interactions.    For  this  purpose  we  introduce  in  the  next  section  a  more 
compact  tetradic  notation.    Furthermore,   it  should  be  noted  already 
at  this  stage  that  we  will  interchange  the  sequence  of  approximations 


8 


with  respect  to  paper  II  by  deriving  the  generalized  unified  theory 
before  making  the  no  quenching  approximation.     This  makes  the  results 
of  the  unified  theory  also  useful  for  situations  where  the  no  quenching 
approximation  cannot  be  made  like,  for  example,  microwave  lines. 
III.     THE  TETRADIC  NOTATION 
The  purpose  of  the  tetradic  notation  which  we  shall  use  is  to 
write  the  product  of  the  U  (t)  operators  in  Eq.  (II.  13)  in  terms  of  a 
single  operator.     To  do  this  we  first  consider  the  product  of  the  matrix 
elements  (a  I  aI  a'  )  and  (  B  I  B!  g  ' )  where  A  and  B  may  be  any  arbitrary 
operator.     This  product  may  be  written  in  terms  of  the  direct  product 
A|^B  according  to 

(alAla')  <e!B|e'>  =  (a^\ABB\a'B')  ,  (IILl) 

where  the  product  states  laB)  =  la)lB)  are  essentially  the  same  as  the 
states  of  Barangers  "doubled  atom" (Baranger ,  1962).    This  direct 
product,  Ag)B,  is  a  simple  form  of  tetradic  operator.    K  one  of  the 
operators  A  or  B  happens    to  be  a  unit  operator  I,  we  may  conveniently 
denote  this  fact  by  means  of  superscripts  t  and  r  according  to 

<aSl  A2;i|a'B')=  <aB  !  A'^  la 'B ')  =  (a\A\a')  (III.  2) 

Bp 

<aBll^-Bla'p')  =  (aSlB^'la'B')  =  <b!bIb')  5     ,    .  (III.  3) 

act 

That  is,   a  superscript  ^  denotes  a  "left"  operator  which  operates  only 
on  the  "left"  subspace  (in  this  case  the  la),   la')  subspace)  and  a 
superscript  r  denotes  a  "right"  operator  which  operates  on  the  "right" 
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subspace.  It  is  thus  clear  that  any  "left"  operator  will  commute  with 
any  "right"  operator: 


[A^,B^]  =  0. 


(III.  4) 


With  this  notation,  the  thermal  average  in  Eq.  (II.  13)  can  now  be 
written  in  the  more  compact  form 


<blu"^{t)lc>  (dIU  (t)|a) 
a  a 

(cd|  TJ^  "\t)  U'^jba) 


av 


(III.  5) 

(clu''(t)|b>  (dIU  (t)la> 


av 


av 


=  (cd| 


U      (t)  (t) 

a  a 


av 


I  ba>  . 


con- 


We  have  chosen  to  write  (blu'^(t)Ic)  as  <  c  I  U  ='T(t)  j  b)  simply  for 
venience  in  the  derivation  given  in  later  sections.     Noting  the  definition 
of  U^(t)  given  in  Eq.  (II.  12),  we  define  operators        (R,x,  v,  t)  and 
V     (R,  X,  v",  t)  so  that 


(R,  ^,  V,  t)  =      exp  I  - 


t 

.     ,         (R,5^,  ^,t')  dt' 


(III.  6) 


U   (R,  5c,  V,  t)  =  O  exp  )         /  V  "^(R,  ^,  ^,  t 


•t 


')  dt'j 


Since  any  "left"  operator  commutes  with  any  "right"  operator, 
have 


we 


U^;^"~(R,  t  V,  t)  U''(R,  X,  V,  t)  =  <^  exp 
a  a 


=^(R,x,i;,t) 


(in.  7) 
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where 


r{R,t,v,t)=  V""  (R,x,  V,  t)  -  V^''(R,x,  V,  t).  (III.  8) 

e  e  e 

We  have  now  succeeded  in  replacing  the  two  U  (t)  operators  by  a  more 
general  tetradic  operator  "^(t)  which  operates  in  both  "left"  and  "right" 
subspaces.    Eq.  (III.  5)  thus  becomes 

r<blU^t)lc)  <d|U  (t)la>l       =  <cdim(t)l  Iba). 

(III.  9) 

It  is  important  to  realize  that  the  tetradic  operator   ^(t)  is  formally 
the  same  as  the  operator  U  (t);  that  is,  it  satisfies  the  same  type  of 

3. 

differential  equation 


ih^^(R,5,v,  t)  =  r^  (R,x,  v,t)^(R,x,v,t)  .  (III.  10) 

This  means  that  all  of  the  line  broadening  formalism  which  has  been 

developed  for  U  (t),  will  be  directly  applicable  to  ^^(t). 

To  make  the  formal  correspondence  more  complete  we  use  the 

operators        ,  h"",  V  (R,  X,  V,  t)  and  V  (R,  x,  v,  t)  to  define  the  tetradics 
o        o       e  e 

K     and  y  (R,  x,  v,  t)  according  to 
o  e 

-K:    =        -  H^^  (III.  11) 

o        o  o 

r  (R,x,  V,  t)  =  v'"(R,x,  V,  t)  -  V^^'^  (R,x,  V,  t)    ,  (III.12) 

e  e  e 
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Since  any  left  operator  commutes  with  any  right  operator  we  have 

(t)  =  exp  ;  itHVh  ^  V^(t)  exp  |-itH^/h>         (III.  13) 
^       o      I     e  I         o  I 

I  it^cyh  I  v^{t)  exp  I  -i«ryh| 


=  exp 


Hence 


(III.  14) 


rjR,  X,  V,  t)  =  exp   I  i^^qA  j  "^J^'  ^)  |"^^o^^| 

which  is  formally  the  same  as  Eq.  (II.  11),    It  is  also  obvious  that  both 

V    and  V    will  be  given  by  a  sum  over  binary  interactions  \    or  7/" 
e  e  11 

just  as  in  Eq.  (II.  4). 


1r  (R,  x,  V,  t)  =  (R,  X  .,  v.,  t)  III.  15 

e  4-^   1         J  J 

J 

r.(R,  x.,v.,t)=  vJ;(R,x.,  v.,t)  -  v'[^''(R,x.,  v.,t)    (in.  16) 

The  formal  similarity  between  the  operators  H  ,  V  (t),  V  (t),  U  (t), 

o       e  e  a 

etc.  and  the  tetradic  operators  K^,  ^g^*)'  ^\(^)>  will  greatly 

simplify  the  treatment  of  the  thermal  average  foi*  the  general  case  of 
upper  and  lower  state  interactions. 
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IV.    THE  GENERALIZED  UNIFIED  THEORY 
Using  the  tetradic  operators  as  defined  in  the  previous  section 
jj  we  have  for  the  correlation  function 

C(t)=    2^  (aldlb)  <cld!d>  e  """"(alp^la) 
abed 

<cd|"F(t)  I  ba)  (IV.  1) 

where  5  (t)  denotes  the  thermal  average  of^*|^(R,  x,  v,  t): 

=  yd^dv  P(x)W(v)^(R,  X,  V,  t)    .  (IV.  2) 

This  tetradic  operator  3  (t)  is  formally  identical  to  the  operator  F(t) 
defined  in  Sec.  (2.  A)  of  paper  III.    It  would  also  be  formally  identical 

y  to  the  F(t)  defined  by  Eq.  (19)  of  paper  II  if  we  would  make  the  no 

I 

quenching  approximation  at  this  point.     To  preserve  generality,  however, 
the  no  quenching  approximation  will  be  deferred  until  a  later  section 
when  we  specify  the  la),  lb),  *  *  *  eigenstates  to  be  H^  eigenstates  for 
hydrogen. 
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Following  the  formalism  developed  in  Sec.  2  of  paper  III,  we 
define  an  operator  ^{H,  x,  v,  t)  by 

g;(R,x,v:t)=  P(x)  W(v)^(R,x,  ^?,t)  (IV.  3) 


so  that 


:3{t)  =   /  d5"dv  3;(R,  x,  v,  t)  (IV.  4) 


(cf.  Eqs.  (11 )  and  (12)  of  paper  III).    From  Eq.  (III.  1  0)  we  see  that 


-  „  .  ih  ^  5  (R,  %  v,t)  =  r        X,  V,  t)  3?(R,  X,  ^,  t)  (IV.  5) 

o  t  e 

which  is  formally  the  same  as  Eq.  (13)  in  paper  III.    We  next  introduce 
a  projection  operator  P  which  is  identically  the  same  as  the  operator  P 
defined  by  Eq.  (14)  of  paper  III  (the  fact  that  P  now  operates  on  tetradics 
does  not  change  its  definition).     That  is,  for  any  function  of  electron 
variables   f(x,  v)  we  have 


Pf(x,  v)  =  P(x)  W(^/)ydx'dv'  f  (x',  v')     ,  (IV.  6) 

This  relation  holds  whether   f   is  a  matrix,  tetradic  or  any  other  type 

of  operator.    With  this  operator  we  can  follow  the  derivation  in  Sec.  (2.  B) 

of  paper  III  replacing  H  ,  V  ,  V    etc.  hyK  ,">/  ,'if    ,  etc.    As  a  result 

o     e      e  o     e  e 

(cf.  Eq.  (27)  in  paper  III)  we  have 
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r^r  3(t)  =  -  h"^/  exp  (itV  A)    y   (t-t')  Q(t-t')r  (0) 
ht  X  o        L  e  e      J  av 


exp  {-it'K  /h)  3(t')  dt'  (IV.  7) 

o 


where 

•t-t' 


Q(R,x,^,  t-t')  =  O  exp  |-^^(l-P)r^{R,x,  V,  t")  dt"|  .        {IV.  8) 

Returning  to  Eqs.  (II.  5)  and  (IV.  1)  we  see  that  the  quantity  of  interest 
is  not  3(t)  but  rather  its  Fourier  transform. 


/•«  .      -iju  t 
(cdlJ(aj)lba)  =/    e"^  e  ( cdl W(t)l  ba>  dt 

Jo 

e  ^^(cdl  exp  (-ittr  /h)'^(t)Iba)  dt 


=J  e"^^  <cdl  3^t)Iba)  dt  (IV.  9) 


where 


g;(t)  =  exp  (-iuc  /h)  3(t)  .  (IV.  10) 

o 


From  Eq.  (IV.  7)  we  see  that 


-|-?(t)  =  -(i^-  /h)'?(t)-h"7  exp  (-i(t-t')ar  /h) 

ht  o'  o' 


I    exp  (-i( 
•'o 

(t-f)  Q(t-t')r  (0)1      ^(t')  dt'.  (IV.  11) 

j  e  Q      I av 


Solving  this  equation  by  Fourier  transforms  gives 


-9(uu)  =  i    I  Auu     -£(Auu     )\  (IV.  12) 


op  op  J 
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where 


X{AUJ     )  =  -  ih 
op 


exp  (itAuu    )   r  (t)Q(t)r  (0) 


op 


av 


dt 


(IV.  13) 


and  Auu 


is  an  operator  defined  by- 


op 


tm     =  ud-k:  /h  =  aj-(H^-H^-'  )/h  . 
op  o  o  o 


(IV.  14) 


With  these  results,  the  line  shape  given  in  Eq.  (II.  5)  becomes 
(cfEq.  (I.  D) 

I((JU,e.)  =  —  Im  (a!dlb)  <c]dld)  (ajp  la) 


(cdl     A«3     -Z^m     )1    ^  I  ba)  . 
[op  op  J       '  ' 


(IV.  15) 


We  next  simplify  jb(Auu^p)  by  means  of  the  impact  approximation 
(see  Sec.  (3.2)  of  paper  II).    Basically  this  approximation  assumes 
that  the  average  collision  is  weak,  that  strong  collisions  do  not  overlap 
in  time  and  that  a  weak  collision  overlapping  a  strong  one  is  negligible 
in  comparison  (weak  collisions  are  those  interactions  for  which  a  low 
order  perturbation  expansion  in  T/"    provides  a  good  approximation  to 

or  Q    ;  for  strong  collisions  the  full  exponential  must  be  retained). 
It  should  be  emphasized  again  that  we  make  a  distinction  between  the 
impact  approximation  and  the  impact  the'ory.      The  latter  contains  the 
impact  approximation  as  well  as  other  approximations  like  the  com- 
pleted collision  assumption  which  will  not  be  made  here.    We  also 
assume  that  the  electron  perturbers  may  be  replaced  by  statistically 
independent  quasi  particles  (e.  g.  ,   shielded  electrons).    In  Sec.  (3) 
and  Appendix  B  of  paper  III,   it  is  shown  that  these  approximations 
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reduce  ^(Atu     )  to 
op 


•t  itAuD 


e       °^  3     (t)  dt  AO)  {IV.  16) 

op 

u 


where 


^\t)=  d5?^dx^^W(v^)  j^<^^(R,5rJ,^?^,t)  -  ij    (IV.  17) 


and 


^^(R,  x^,  v^,t)  =  (yexp  Z^j(R,  x^,  v^,t')dt'|(IV.  18) 

and  n    denotes  the  electron  density, 
e 

Equations  {IV.  15)  through  {IV.  18)  give  the  line  profile  of  the 

generalized  unified  theory.     To  obtain  the  impact  theory  we  simply 

replace  i!{Auu     )  by  Z{0)  and  as  discussed  in  Sec.  4  of  paper  III,  we  have 
op 

the  familiar  result  (cf.  Eq.  (44)  of  Baranger,  196Z). 

£(0)  =  i y [S^'"  S^^  -  1)  dv  (IV.  19) 


where  S^^  denotes  an  S-matrix  for  a  binary  (completed)  collision  and 
^ dv  denotes  the  integral  over  collision  variables,  as  defined  in  the 
Appendix  of  paper  II. 

00 

ydv  -  n^  ^  dv  vfCv) dp  2rrp J' dfl  (IV.  20) 

In  comparing  Eq.  (IV.  19)  with  Barangers  result  it  is  important  to  note 
that  Barangers  operators  S.  and  operate  only  on  "initial"  and  "final" 
states  respectively,  whereas  our  operators         and  S^^    operate  on  all 
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possible  H    eigenstates.    This  difference  occurs  because  we  have  not 
o 

made  the  no  quenching  assumption  yet. 

The  other  limit  of  the  one  electron  theory  is  obtained  by  making 


a  wing  expansion  of  the  unified  theory;  that  is,  the  operator 

[Am     -I  (Acu  )]' 
op  op 

order  this  gives 


[A(«  (Acu     )1   ^    is  expanded  in  powers  of  [jb(Atw     )/A(A)     1.    To  lowest 

op  op  op  op" 


[m    -ziisw    )'\~^  = + zim    )-r-  +  

op  op  A.UJ        Auu  op  A^ 

op        op  op 

■I  rt  ^^^^      n ) 

=  ^      -  i/    e  V)dt+  •••    .  (IV.21) 

op  •/o 

The  first  term,   ^/A^u^p,   gives  a  delta  function  when  one  takes  the 

imaginary  part  required  by  Eq.  (IV.  15).     To  get  this  delta  function  we 

approximate  radiation  damping  effects  by  using  (a^d      +  ie)  in  place  of 

op 

AOJ      (see  Sec.  (3.  A)  of  Smith  and  Hooper,  1967);  the  imaginary  part  of 
op 

1  /Muu      +  ie )  is  just  -nbiAiU     )  when  e->0.    When  this  delta  function  term 
op  op 

is  averaged  over  ion  fields  according  to  Eq.  (II.  1)  it  will  produce  the 
line  broadening  due  to  the  static  ions  alone  (see  Sec.   5  of  paper  II). 
The  influence  of  the  electrons  as  well  as  electron-ion  coupling  is 
contained  in  the  second  term  of  Eq.  (IV.21).    Hence  one  is  interested 
in  the  matrix  elements  of  the  Fourier  transform  of  3^^^(t),  which  is 
also  the  quantity  of  interest  in  the  unified  theory  (see  Eq.  (IV.  16)).  The 
primary  difference  between  calculations  made  by  the  unified  and  one- 
electron  theories  is  therefore  the  matrix  inversion  of  f  Auu     -X(Auu  )1 

op  op 

which  is  required  by  the  unified  theory  but  not  by  the  one  electron 

-(1 ) 

theory.    Since  the  matrix  elements  of  the  Fourier  transform  of  3"  (t) 
play  such  a  central  role  in  any  classical  path  theory  (including  the  im- 
pact theory),  the  evaluation  of  these  matrix  elements  for  hydrogen  will 
be  discussed  in  detail  in  the  following  sections. 

18 


V.    THE  NO-QUENCHING  APPROXIMATION  FOR  HYDROGEN 

In  the  preceding  section  we  have  derived  the  thermal  average 
3f{t)  and  its  Fourier  transform  .9  (uo )  for  the  general  case  of  upper  and 
lower  state  interaction.    In  order  to  evaluate  I((jo,S.)  in  Eq.  (IV.  15) 
we  have  to  consider  the  complete  trace  over  all  H  eigenstates 
la),   lb),  •••  .    However,  in  looking  at  the  Eqs.  (IV.  9),  (IV.  18)  and 
(III.  14)  one  realizes  that  due  to  the  exponential  factors  only  a  few  of  all 
the  possible  matrix  elements  will  contribute  significantly  to  the  final 
line  profile  at  a  particular  frequency  uu.     That  is,  we  can  neglect  those 
matrix  elements  for  which  the  argument  of  the  exponential  factor  is 
so  large  that  it  gives  rise  to  rapid  oscillations  within  the  range  of  the 
time  integral.    Hence,  if  one  treats  well  isolated  lines,  only  those 
matrix  elements  of  Uj(t)  between  either  "initial"  or  "final"  states 
have  to  be  considered.    We  may  therefore  state  the  no-quenching 
approximation  as 

^j(t)=  V^^\t)\j\{t)  (V.l) 

where        now  no  longer  operates  on  the  complete  "left"  or  "right" 
subspace,  but  only  on  "initial"  or  "final"  states  (  see  also  Sec.  2.2 
and  7.  2  of  paper  II). 

Further  approximations  cannot  easily  be  generalized  and  depend  on 
the  particular  problem  investigated.    We  now  apply  our  general  results 
to  the  problem  of  hydrogen.    In  this  case  the  no-quenching  assumption 
states  that  we  need  to  consider  only  those  matrix  elements  of  U^^t) 
and  '\^^(R,  ^jj'^j'  ^)  which  are  diagonal  in  the  principal  quantum  number 
n.     As  shown  already  in  paper  II  this  is  a  good  approximation  as  long 
as  the  lines  investigated  are  well  separated.    For  calculating  the  line 
wings  it  is  furthermore  required  that  there  is  no  appreciable  overlap 
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with  wings  of  adjacent  lines  in  the  region  of  interest.     The  same  is 

true  also  in  any  reliable  measurement  of  line  wings. 

To  show  this  we  can  proceed  as  in  Sees.  2.  2  and  7.  2  of  paper 

II  with  the  difference  that  now  we  are  dealing  with  the  operator 

H    =  H    +  eZ£.    rather  than  just  H  .    Since  H    does  not  commute  with 
o        a  1  a  a 

Z  we  introduce  a  projection  operator  P    (see  Sec.  2.  2  of  paper  II) 

n 

which  picks  out  the  part  of  an  operator  which  is  diagonal  in  n.  Using 
this  operator  we  split  H    into  a  part  which  is  diagonal  in  n 


H      =  H    +  eP  Ze.  (V.  2) 

on        a  n  1 


and  a  part  which  is  not  diagonal  in  n 


H  .  =  e{l-P  )  Ze..  (V.  3) 

off  n  1 


H    now  commutes  with  P   Z  because  both  operators  are  diagonal  in 
a  n 

parabolic  Coordinates.    We  therefore  specify  their  eigenstates  com- 
pletely by  the  principal  quantum  number  n,  the  magnetic  quantum 
number  m  and  the  quantum  number  q  which  is  defined  to  be 


q^n^-n^;  (V.  4) 


n  and  n  are  the  usual  parabolic  quantum  numbers  which  obey  the 
relation 


n  =  n    +  n    +  I  m|  +  1   .  (V.  5) 
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Knowing  the  solution  of  the  eigenvalue  problem 


H       nqm)  =  E         I  nqm)  (V.  6) 

on'  nqnn' 


with 


E         =  E    +  eZ    p,.  (V.  7) 

nqm        n  nq  i 


we  see  from  a  second  order  perturbation  approach  (cf.  Chap.  1  6  of 
Merzbacher,   1961)  that  the  energy  correction 

2 


(2)       ^  j<nqm!H^^^ln'q'm')| 


I  ^^nqm=2.,^-^  -E   .    .     .  (^-^^ 

rFn*  nqm  nqm 

can  always  be  neglected  as  long  as  the  ion  fields  do  not  become  too 
large.    This  is  again  equivalent  to  stating  that  the  lines  have  to  be  well 
separated. 


As  a  result  one  is  left  with  the  eigenvalues  E  ,   E  ,   

a  b 

of  the  Hamiltonian  H      whose  eigenstates  la),   lb),  •  •  •  I  d)  are  the 

on 

parabolic  states  I  nqm) .    This  allows  us  to  rewrite  the  autocorrelation 
function  C(t)  in  Eq.  (IV.  1)  for  hydrogen  in  the  form: 

C(t,r..)  (nq^m^l  dl  n'q^m^  )  <  n'q^m^l  dl  nq^mj^  ) 

exp  f-  —  [e  -E   /+e{Z       -2   ,    ,   )e.lt|(nqm  !p  Inqm) 
^  y     h    L  n     n'  nq^     ^^b  )        ^    a     a       a  a' 

^n'q'  m'  ;  nq^m^!?{t)!n'q'm'   ;  nq  m   )  (V.  9) 

bb        bb  aa  aa 

where  quantum  numbers  which  refer  to  the  lower  state  are  distinguished 
from  the  upper  state  quantum  numbers  by  a  prime. 
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The  matrix  elements  of  P  Z  are  given  by  (see  Bethe  and 


n 


Salpeter,  1957) 


(nqmlzlnqm)  =  Z      =   —   nqa  (V.  10) 

■  nq       2  o 


2  2 

with  a^  =  h  /(me  )  being  the  Bohr  radius.    As  a  further  definition  the 

ion  field  6.  will  be  normalized  to  the  Holtsmark  field  strength  6 
1  ^  o 


P..  =  B-e  (V.  11) 

1  o 


where 


■      (  4tjY 


2 

l-T")   e  n^    .  (V.  12) 


This  yields 


-1-{Z       -Z   ,  =  AUD.(n,q  ,n',q')  •  R  {V.13) 

h   ^    nq^     n'q^     lib  b 


with 


/        -  - 
Aa).(n,  q^,  n',  q^)  =  .  -\   (nq^-n^q^)  ^  n^  . 

(V.  14) 


Am.  is  now  the  frequency  shift  of  a  particular  Stark  component  char- 
acterized by  the  quantum  numbers  ^>  "^^j  ^'  ^^'^  to  the  Holtsmark 
field  strengths  .    Introducing  the  frequency  shift  t\\Xi  =  ,  where 


nn 


the  frequency  of  the  unperturbed  line  uu      ,  is  given  by 


nn 


(D     ,  =  (E  -E   ,)/h  (V.  15) 

nn  n  n 
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the  line  profile  I(Auu,  p)  can  be  written  in  the  form 


I(AU3,B)  =  —  y  <nq  m  I  dl  n'q'm')  <n'q'mM  d|  nq  m  ) 
TT  ^— ^       a    a  a    a  bbbb 


(V.  16) 


(nq  m  Ip    I  nq  m  )  ('n'q'm' ;  nq^m,  |  j(tju  )|  n'q'm' ;  nq  m  ) 
a    a     a       aa  bb       bb'  aa  aa 


where 


(n'q'm'';  nq  m  L9(uu)  In'q'm';  nq  m  )  = 

bbbb  aaaa  ^^^^^^ 


X 


dt  exp  <  i(Acu -Auu  .B  )t  >  (n'q/m/ ;  nq,        I  5(t)l  n'q'm' ;  nq  m  )  . 
^  i|  bbbb  aa  aa 


Performing  the  ion  field  average  according  to  Eq.  (II.  1)  will  then  give 
us  the  desired  line  profile  once  we  know  the  thermal  average  ^(t). 
VI.    THE  THERMAL  AVERAGE  ^^^V)  FOR  HYDROGEN 

In  Sec.  ,IV.  we  saw  that  the  crucial  problem  in  any  classical 


path  theory  of  line  broadening  is  the  evaluation  of  the  matrix  elements 
of  3^^^(t).  With  the  no-quenching  approximation  for  hydrogen  a  typical 
matrix  element  in  parabolic  states  is  given  by 


(n'q'm';  nq,  m_  I  3?^^  ^(t)]  n'q'm' ;  nq  m  >  = 

bb       bb  aa       a    a  ^  ^ 

^/^^l^^l^^^l^     (n'q^m;  ;  nq^m^l^2/|^(t).l  ]  n'q^m^;  nq^m^> 

To  simplify  the  evaluation  we  transform  to  the  natural  collision  vari- 
ables p,  V  and  t    which  denote  the  impact  parameter  electron  velocity 
o  ' 

and  some  reference  time  of  the  collision  (see  the  appendix  of  paper  II). 
The  orientation  of  the  collision  axes  with  respect  to  the  radius  vector 
R  of  the  orbital  electron  is  specified  by  the  three  Euler  angles 
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represented   by  Q.    Furthermore  we  assume  a  spherically  symmetric 
distribution  of  perturbing  electrons;  this  is  a  good  approximation  as 
long  as  the  impact  parameters  are  not  too  small.    The  velocity  dis- 
tribution function  W(v)  is  related  to  the  Maxwell  distribution  function 
f{v)  by 


f{v)  =  4ttv^W(v).  (VI.  2) 


With  the  preceding  definitions  Eq.  {VI.  1)  can  be  rewritten  as 


(n'q'm'   ;  nq  m,  j'S^^  V)l  n'q'm';  nq  m  )  = 
b      b       b    b  a    a       a  a 


{VI.  3) 


^Jdnjdv  vf(vydp  p^dt^  ^^'^^{^^^J\}^ithi\n'<^^m^^ 

Next  we  have  to  know  the  matrix  elements  of  the  time  development 
operator  Ulj^{t)  defined  by  Eq.  (IV.  18).    This  requires  the  matrix 
elements  of  the  interaction  potential  ^^(t).    In  order  to  save  some 
writing  we  consider  for  the  moment  only  Uj^(t)  and  Vj^(t)  which  after 
making  the  no-quenching  assumption  may  be  the  "initial"  or  "final" 
part  of  the  corresponding  tetradic  operators  (see  Eq.  (V.  1)).  A 
typical  matrix  element  of  Vj^(t)  is  given  by 

<nq  m  I  V,  (t)l  nq  ,m  ,)  =  expJi—  (Z  -Z  )e.tWnq  m  j  V,  (t)j  nq  ,m  , 
^   ^c    c'    r    '   ^d    d'^  h        nq       nq  /   i  I  ^   ^c    c'    1      '    ^d  d 

^  c  d         )  ^^j^^j 


With  the  no-quenching  assumption  the  unperturbed  energy  eigenvalues 

E    have  cancelled.    At  this  stage  we  now  make  another  simplification 
n 

by  dropping  the  exponential  in  the  latter  equation;   this  has  been  done 
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jin  all  previous  Stark  broadening  calculations  but  it  is  rarely  stated 
explicitly.    This  will  be  a  good  approximation  in  the  line  wings  where 
the  times  of  interest  1/auu  are  small  and  Auu  is  much  larger  than  the 
average  ion  field  splitting.    In  the  line  center,  however,  the  argument 
of  the  exponential  can  easily  be  on  the  order  of  unity  or  larger  in 
which  case   V^(t)  effectively  vanishes  due  to  rapid  oscillations  of  the 
exponentials.     This  effect  was  first  noted  by  Van  Regemorter,  1 964  ,  who 
shows  that  this  effectively  introduces  another  cut  off  which  may  easily 
be  smaller  than  the  usual  Debye  or  Lewis  cutoffs  .    This  additional 
cutoff  has  been  included  in  recent  calculations(  Kepple  and  Griem,  1968). 
However,     as  discussed  in  Sec.  XII  it  turns  out  that  its  influence  on  the 
final  line  profile  is  in  most  cases  negligible. 

Neglecting  the  ion  field  exponentials  in  Eq,  (VI.  4).  the  time 
development  operator  U^^  is  now  given  by 

where  the  time  ordering  is  still  required  because  P  Vj^(t)  need  not 
commute  with  P  V^(t').    In  paper  II  it  was  shown  that  this  time  ordering 
is  negligible  for  weak  collisions  (to  second  order)  as  well  as  quasi- 
static  collisions  (i.  e.  ,   in  the  distant  line  wings).     Time  ordering  is  not 
negligible  for  strong  collisions;  however,  when  the  thermal  average  is 
performed,  the  errors  due  to  neglecting  time  ordering  are  expected  to 
be  small.     The  reason  for  this  is  that  the  time  development  operator 

still  retains  its  unitarity  (cf.  Sec.  8  of  paper  11). 
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VII.     THE  MULTIPOLE  EXPANSION  OF  THE  CLASSICAL 
INTERACTION  POTENTIAL 
Before  evaluating  the  thermal  average  5^^\t)  we  briefly  consider 
the  classical  interaction  potential  V^(t)  due  to  a  single  electron.    If  the 
perturber  does  not  "penetrate"  the  radiator,  Vj^{t)  is  given  by  the  well 
known  multipole  expansion 

^  [r(t)1 

where  I  R|  is  the  distance  of  the  orbital  electron  fronn  the  nucleus,  r(t) 


00  ^  ]^ 

V,(t)  =        T;  — •  Pk  [cos  e(t)l      (VII.  1) 


is  the  instantaneous  distance  of  the  perturbing  electron,  the  P  are 
Legendre  polynomials  and  0(t)  is  the  instantaneous  angle  between  R 
and  r(t). 

In  most  cases  it  is  sufficient  to  consider  only  the  dipole  (k=l) 

term.    However,  to  account  for  some  asymmetries  of  a  line,   it  may  be 

necessary  to  keep  some  of  the  higher  multipole  terms  as  well.    In  any 

case,  one  can  show  that  this  multipole  expansion  is  terminated  after 

some  finite  number  of  terms  due  to  symmetries  of  the  radiator. 

->  i> 

To  show  this  we  specify  the  angular  positions  of  R  and  r(t)  by 
9.,,  cp.,  and  0  ,  cp    respectively  and  we  apply  the  spherical  harmonic 
addition  theorem  (Eq.  (4.6.7)  of  Edmonds,  I960) 

+k 

P^(cos0)=^   (-1)^       (0^,  cpj)  •  0^02, cp^)  (vn.2) 
p=-k  ^  ^ 

where 


cos  0  =  cos  0    COS  0    +  sin  0     sin  0     cos  (cp  -cp_).  (VII.  3) 

X  Ld  X  ui  L  iL 
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We  may  simplify  the  mathematics  without  loosing  generality  by  choosing 
a  coordinate  system  in  which  cp    =  0.    Using  the  relation 


(VIL  4) 


one  then  obtains 


p  fcos  e)=     (e,,cp,)c^  (9,)  + 

k  o     1     1     o  Z 


p=i 


which  gives  for  the  interaction  potential 


\tl^ 


k=l  [r(t)] 


C    .  P,  (cos  e-(t))  + 
ok  d 


The  dipole  case  {k=l)  gives  the  well  known  result 


r  (t) 


cos  9,(t)+^jc^j-C;}si 


(VII.  5) 


(VII.  6) 


sin  Q^{t} 


Z  •  cos  Q^{t)  +  X     sin  B^{t) 


] 


(VII.  7) 


The  y- component  vanished  because  cp  =  0.    Similarily  one  can  write 
down  the  higher  order  multipole  terms.    The  necessary  matrix  elements 
of        are  given  by 


P 


/  /-t'k^\/f,'k/f\ 


27 


From  the  last  3j- symbol  we  see  that  these  matrix  elements  exist  only 

if       k,  I,'  satisfy  the  triangle  condition  and  their  sum  is  an  even  integer. 

Therefore  it  turns  out  that,  within  the  no-quenching  assumption  where 

one  needs  the  matrix  elements  of  P  V,  (t),  only  a  finite  number  of 

n  1 

multipole  terms  exist.     The  summation  index  k  in  equation  (VII.  6)  has 
to  obey  the  condition 

1  <  k  <  2  (n-1)   .  (VIL  9) 

As  an  example  we  see  that  a  calculation  of  the  upper  state  interaction 

of  Lyman  a  requires  only  the  dipole  and  quadrupole  terms.  This 

condition  also  illustrates  the  well  known  fact  that  there  is  no  ground 

state  interaction  for  the  Lyman  series. 

VIIL     THE  SPHERICAL  AVERAGE  OF  THE  TIME 

DEVELOPMENT  OPERATOR  Ul  (t) 

_(1) 

In  our  evaluation  of  the  thermal  average      5     {t)jdefined  in 

Eq.  (VI.  3)jWe  first  perform  the  spherical  average  represented  by  the 

integral  over  the  Euler  angles  Q,  because  it  greatly  simplifies  the 

remaining  integrals  over  t  ,  p  and  v.     This  is  due  to  the  spherical 

o 

symmetry  of  the  time  development  operator  Uj(t)  defined  in  Eq.  (VI.  6). 
It  should  be  noted  that  this  symmetry  was  achieved  by  dropping  the  ion 
field  exponentials  in  Eq.  (VI.  4),  thus  replacing  Vj^(t)  by  V^(t).  We 
will  perform  this  average  by  means  of  a  rotation  technique  used  by 
J.  Cooper,  1967,  andBarangei,   1958,  for  S-matrices.    Although  we 

are  working  with  the  more  general  time  development  operators 
Uj^(t)  or ']j^j(t),  the  rotation  technique  is  the  same. 

In  terms  of  the  collision  variables  p  ,  v,  t    and  Q,  the  dipole 

o 

interaction  between  the  radiator  and  a  perturber  is  given  by 

3 

^  (VIIL  1 ) 


V^(t)  =  e  R 


p  +  V  (t+t  ) 
o 


.2  2,  ,2 
/   p  +v    t+t  ) 

'  o 
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(see  the  appendix      of  paper  II).     The  three  Euler  angles  denoted  by 

describe  the  orientation  of  the  collision  frame  relative  to  the  atomic 
frame.    It  is  therefore  convenient  to  perform  a  rotation  of  the  atomic 
axis,  through  the  angles  Q  in  such  a  way  that  R  points  in  the  same  dir- 
ection as  ^  and  the  x  axis  of  the  rotated  atomic  frame  points  in  the 
same  direction  as  "v.    In  this  rotated  frame,  the  interaction  potential 
takes  the  form 


V  (t)  =  e^ 
c 


Zp  +  Xv  (t+t^)j  /  j^p^+v^{t+t^)^  j  ^  (VIII.  2) 


I  This  rotation  transforms  the  time  development  operator        into  a  new 

I  operator  U,    ,  where  U    and  U,     are  related  to  one  another  by 
^  Ic  1  Ic 

\J    =    &  ~^(Q)  m)  (VIII.  3) 

1  1  c 

I  where  &  (Q)  is  a  rotation  operator  (see  Chap.  4  of  Edmonds,  I960).  The 

j 

i  time  development  operator  in  the  rotated  frame,  U^^,   is  given  by 


U,     =  exp 
1  c 


!-  -  IP  V  (f)  dt'  >    .  (VIII.  4) 


To  make  the  form  of  U,     more  explicit,  we  perform  the  integral  over 

1  c 

t'  and  we  obtain 


J     h     pv  1^  n 


U,     =  exp  <-  —  —  P  ZA(t,  t  ,p,  v)-P  XB(t,  t  ,p,  v)J    (VIII.  5) 
Ic  JhDv       n  o  n  o  * 


(v/p)(t+t)  (vt/p) 
A(t,t^,p,v)=  °  -  .  (VIII.  6) 


V  l+(v/p  )^(t^+t)^  ijlHvtJpy 
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and 


B(t,  t^,  p,v)=    -  —  .      (VIII.  7) 


^fl7(v/p  )^(t^+t)^       ^1  +  (vtyp)^ 


Substituting  Eq.  (VIII.  3)  into  Eq.  (VI.  3)  we  see  that  the  integral  over 
Q  in  Eq.  (VI.  3)  involves  only  the  matrix  elements  of  four  rotation 
operators.    Since  it  is  convenient  to  use  spherical  states  I  nJf.m)  when 
taking  matrix  elements  of  j9  (Q)  ,  we  make  use  of  the  unitary  trans- 
formation from  parabolic  to  spherical  states  discussed  by  Hughes,  1967. 


I  nqm)  =        I  nf.m '  )  (  n^,m '  j  nqm) 


/    ,      /I         \       c  /   ,  .^(1+m-q-n)  /      .  , 

(rLf.m' I  nqm)  =  6        ,(-1)  J  ^1+^ 

mm  % 


2  2 


using  3j- symbols  and  the  definitions  in  the  Eqs.  (V.  4)  and  (V.  5).  (An 

error  in  the  phase  factor  has  been  pointed  out  by  H.  Pfennig,  private 
communication).    Noti-ng  that  S5(0)  is  diagonal  in  the  angular  momentum  Jf,, 
the  Q  integral  in  Eq.  (VI.  3)  may  now  be  written 


y(n'q'm';  nq  m   |  Ul  (t)|  n'q 'm ' ;  nq  m  >  dQ  = 
bb       bb       l  aa  aa 


(^')-l 


E/dd(n'q'mMn'^'m')   J5  "  ^'       (n'^, 'm' I  U^"'"  I  n'^ 'm') 
J    \        aa        aa  aclc  bd 

"  mm 

^    ^  (VIII.  9) 


no 


^      ,     ,  (n'^'m'ln'q'm') 
m'm'   ^     ^b    b         b  b 
d  b 

m 


(nq^m^ln^^n:^) 

1  m,  bd  lcac  aaaai 
h    d  mm  J 
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where  the  summation  7^  denotes  sums  over  ^  ,  ^  '  9^     rn     m'  m 

^  a    a      b '   b '      c'     c'  d 

and  m'.     The  Q  dependence  of  the  integrand  in  Eq.  (VIII.  9)  is  contained 
d 

entirely  in  the  four  rotation  operators  ^  (Q).    Using  Eqs.  (4.  3.2)  and 
(4,  6. 1  )of  Edmonds,  1 960,  we  obtain  the  identity 


/ 


a                 b           b  a  _ 

&  dQ  = 

mm         m '  m      mm  m  m 

ac          db       bd  ca 


9     --^  m'-m'  +  m  -m,  +  M-M' 

S  (-1)    ^      "        ^      ^  (2L+1) 

L,  M,  M' 


(VIII.  10) 


aa  Uaa         i/bb  Vbb 


•  m '    m      My  \m '  m      M /  Vm '    m  ,    M y\m'    m,     M , 
aa  dd  ^bb 


Hence  Eq.  (VIII.  9)  becomes 


J^'^'^^h^W  ''%'^bl'^l^^^'''''^I'^a'  '''^a'^a^  " 


E(n'q'm' 1  n'-f/m')  (n'^/m/ I  n'q/m/')  (nq^  m,  !  nP,,  m,  )  <frLf,  m  |  nq  m  ) 
a    a        a    a  b    b        b    b        b    b      b    b        a    a       a  a 

^   ,  m'-m'  +  m.-m, +M-M' 

8^2  (.1)     =      a        d  b 

(VIILU) 


•m'm  My  Vm'  m   M/  \-m',  m^MyVm'    m,  M, 
cc  aa  dd  bb 


(n'  l'  m'  ;  nf ,  m  ,  !"7n  ,    (t)l  n';/m';  n;   m  ) 
^b     d        b    d     UVlc  a    c       a  c 
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This  result  is  spherically  symmetric;  that  is,  any  further  rotations  of 

the  atomic  coordinate  system  leave  this  expression  unchanged.  One 

may  verify  this  rotational  invariance  by  rotating   U      through  some 

-1 

arbitrary  angle  Q '  so  that  U,     =     ^     (Q')U'      j&(Q').     Taking  matrix 

1  c  1  c 

elements  of  the  new  rotation  operators  and  making  use  of  the  ortho- 
gonality properties  of  the  3j- symbols  one  sees  that  the  right  hand  side 
of  Eq.  (VIII.  11)  did  not  change.    Since  we  are  free  to  perform  further 

rotations  on  U,     without  altering  Eq.  (VIII.  11),  it  is  convenient  to 
1  c 

rotate  the  X-Y  plane  through  an  angle  e  =  arctg  (B/A)  where  A  and  B 


are  given  by  Eqs.  (VIII.  6)  and  (VIII.  7).  This  rotation  transforms  U 
into  an  operator  U^^  given  by 


=  exp 


—  e  P  Zg(t,  t  ,p,v) 
h  n  o 


Ic 


(VIII.  12) 


where 


_L  ^ 

pv 


g(t,  t^,  p,v)  =         yA    +  B- 


pv 


1  +  (v/p)  t    (t  +  t) 
o  o 


(VIII.  13) 


1  -  — 

y  1  +  (v/p)^(t^  +  t)^  ^1  +  (vtyp)^ 


The  operator  U^^  has  the  important  property  that  it  is  diagonal  in 
parabolic  states  (because  it  contains  only  P  Z).    Hence  a  typical 


n 


matrix  element  of  U^^    is  given  by 


(nqml  U^(t)|  nqm)  =  exp    |-i  "|"  ^q  ^  g  (t,  t^,  p,  v)| 


(VIII.  14) 


We  also  realize  that  one  and  the  same  rotation  through  the  angle  e=  arctg 
(B/A)  diagonalizes    simultaneously  both  time  development  operators 
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acting   on  initial  and  final  states  respectively.    As  a  result  a  typical 
matrix  element  of  the  corresponding  tetradic  operator  iJf^(t)  is  given  by 

<n'q'm';nqna|'Ul^(t)|  n'q'm';nqm)  =  exp  | -i -|- (nq-n'q '  )^g't,  t^p  ,  v)| 


(VIII.  15) 


Substituting  this  identity  into  Eq.  (VIIL  11)  the  spherical  average  of  the 
time  development  operator  "U[j^(t)  finally  becomes 

f(n'^'m';  nq  m J  lA  (t)|  n'q 'm'  ;  nq  m  )  dO  = 
/        bb       bb'l  aa        a  a 

X^<n'q'm'l  nV.'m  <n'^ 'm'  I  n'q' m'  >  (n'q'm' I  n'^/m')  (n'^/m/ |  n'q/m/) 
X  y       aa        aa  ac         cc  cc'       be  bb  bb 

(nq^mj  n^^m^Xn^^mJ  nq^m^Xnq^mJ  n^^m^>  (n^^m^I  nq^m^) 


Srr^-l)      ^      Vl+1)(    ^  ^         Y    ^    ^          irb    'b^Y  ^  ^ 

Um^  m^    Myy-m^  m^  MyVm^  m^  M7Vm^  M 

exp  {  -i-|-(nq  -n'qO  —  g  (t,  t  ,p,  v)[  (VIII.  16) 

'        Z        c         cm  o  I 


where  the  unitary  transformations  are  given  by  Eq.  (VIII.  8).  This 
result  greatly  simplifies  if  there  is  no  lower  state  interaction  (e.  g.  , 
Lyman  lines),  in  which  case  one  obtains 


/ 


(nq  m  IU(t)lnq  m  >  dO  =\      (na  m  I  n^,  m  >(n^,  m  !  nq  m  >  5 

D    b     1         a    a  /  J      ^b    b      a    a        a    a      a    a^  m  m 

a  b 

pn^^mjnq^m^)j      f^-^^  exp    j-i^  nq^  ^  g(t.t^,p,v)j 


(VIII.  17) 
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This  simplified  relation  may  also  be  used  for  the  higher  series  mem- 
bers of  the  Balmer^Paschen  etc.  series  where  lower  state  interactions 
contribute  only  a  negligible  amount  of  broadening  to  the  final  line 
profile. 

IX.    EVALUATION  OF  THE  THERMAL  AVERAGE 
"i^^^(t)  FOR  HYDROGEN 
Having  performed  the  spherical  average  over  the  Euler  angles 
Q  we  can  rewrite  Eq.  (VI.  3)  in  the  form 

^'''^b"^b'  "^S^h^  ^     (t)|n'q^m^;  m^rn^)  =  (-1)       ^        ^  (ZL+l) 

<n'q'm'|n'^'m'><n';/m'|n'q'm')<n'q'm'|n'^'m'><n'^'m/In'q/m 
a    a        aa  ac        cc  cc        be  bb  b 

(nq  m  \nl  m  )(rLf,  m  !  nq  m  )(nq  m   I  n^,  m  )  (nt,  m,  I  nq,  m,  ) 
a    a       a    a        a    c       c    c         c    c       b    c        b    b      b  b 

-^c   ^c  k  "^a^/V^c  ^c^/V^b^b  m/^^*'''''^^'''''^^^ 

(IX.  1) 


where 

F(t,  n,  q  ,  n',  q')  =  Znn    /  dv  vf(v) /dp  p /dt      $(t,  t,p,v)       (IX.  2) 
c  c  ej  J       J  ° 

and 

$  (t,t  ,p,v)=  exp  {-i-^  (nq^-n'qM  —  g  (t,  t  ,p,v)\-l.  (IX. 3) 


j-i^(nq^-n'q^)^g(t,t^,p,v)j-: 


Thus,  the  problem  is  now  reduced  to  evaluating  F(t),  which  will  be 
done  in  this  section.    It  is  interesting  to  note  the  similarity  between 
Eqs.  (Vm.   13)  and  (IX.  2)  and  the  i]/ -function  of  Anderson  and  Talman, 
1955,  which  is  the  crucial  function  in  their  classical  adiabatic  theory. 
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We  first  realize  that  due  to  the  symraetry  of  the  line  profile  we 

only  have  to  evaluate  the  real  part  of    $  (t,  t^  ,  p ,  v);  that  is,  for  every 

positive  value  of  (nq  -n'q')  there  will  be  the  corresponding  negative 
^  c  c 

value.    Hence  we  are  left  with 

Mt.t^.p,  v)=  cos  I  ■|-(nq^-n'qM  ^        t^,  p  ,  v)  |  -1 .  (IX. 4) 

In  performing  the  integrals  over  p  and  t^  in  Eq.  (IX.  2)  we  account  for 
shielding  by  setting  the  interaction  potential  V  and  hence  also    $  equal 
to  zero  whenever  the   distance  of  the  perturbing  electron  is  larger  than 
the  Debye  length  D.    We  also  introduce  a  strong  collision  cutoff  p  ^^^^ 

In  principle  we  can  let  the  impact  parameter  go  to  zero  because  the 

functions  $  and  F(t)  do  not  diverge  for  small  impact  parameters  as  they 

do  in  some  second  order  theories.    However,  for  numerical  purposes 

this  would  result  in  very  large  computer  times  due  to  the  growing 

fluctuations  in  the  integral.    For  this  reason  we  will  choose  p     .    to  be 

mm 

small  enough  so  that  when  we  are  interested  in  large  frequency  per- 
turbations Auo  where  perturbers  at  small  impact  parameters  are  quasi- 
static,  the  rest  of  the  integral  from  0  to  p     .    may  then  be  replaced  by 

mm 

the  static  limit.    In  the  dipole  approximation  this  gives  rise  to  the 

-5/2 

well  known  Holtsmark  AX         -wing  (see  also  Sec.  X).  According 
to  the  validity  conditions  of  the  classical  path  theories  (see  Paper  I)  the 

minimum  impact  parameter  p     .    will  be  of  the  order  of 

mm 


2 

p^  =     +  n  a^  (IX.  5) 


where ^  is  the  De  Broglie  wavelength. 
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We  now  concentrate  our  attention  on  the  integral 


/ 


G(t,  p,  v)  =  /  dt    $  (t,  t  ,  p,  v). 
■      o  o 


(IX.  6) 


For  convenience  we  consider  the  collision  sphere  as  shown  in  Fig.  1. 

The  perturbing  electron  moves  along  the  classical  straight  line 

trajectory  L  and  we  are  interested  in  the  interaction  from  some  time 

t  to  some  time  t  +  t.  Due  to  the  Debye  cutoff  the  t  -  integral  extends 
o  o  ^  o  ^ 

from  -T  to  +  T  where 


(IX.  7) 


and  the  interaction  potential    vanishes  if  the  electron  is  outside  the 
sphere  of  radius  D.     The  corresponding  time  integration  limit  t  due  to 


the  strong  collision  cutoff  p     .    is  given  by 
^  ^mm      ^  ' 


1  rr 


(IX.  8) 


Based  on  this  model  of  the  collision  sphere  we  split  the  integral  G 
into  two  parts 


G(t,  p,v)=    U(p>p   )  -G  (t,  p,v)+  U(p  >p)  -Gjt,  p,v) 
o       a  o  b 


(IX.  9) 


where  the  step  function   U   is  defined  to  be 

r 

U(a>b) = 


1       if  a^b 


0       if  a<b 


(IX.  10) 
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In  order  to  evaluate  G  (t,  p,  v)  we  have  to  distinguish  the  following 

3, 

four  cases  depending  on  whether  the  initial  and  final  times  of  inter- 
action are  inside  or  outside  the  sphere. 


Case  I:-T<t;t  +t<T 
o  o 


p.  V 


1  +— 2:  t_  (t  +  t) 


o  o 


14 


(IX.  11a) 


This  is  the  same  general  expression  as  given  in  Eq.  (VIII.  13). 


Case  II:     -T  <  t  ;  T  <  t    +  t 
o  o 


D  p 


P 


(IX. lib) 


Case  III:   t   <-T;   t  +t<T 

o  Op. 


vr 


D  p 


1  — — ^ 


(IX. 11c) 


Case  IV:     t   <  -  T;  T  <  t    +  t 
o  o 


g4(t,  t^,  p,v)  = 


pv 


4 


(IX.  lid) 


D 
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After  defining 


(t,  t  ,  p,  v)  =  cos  I  ^  (nq  -n'q')  g  (t,  t  p,  v)  S 
ko  Id         c         cm.  °k        o  1 


the  integral  G    is  given  by 

3. 


-1 

(IX.  12) 


T-t  T  -T 


G    (t,  p,  v)  =  U(2T>t)|   /$,dt    +/$_dt    +/  $_  dt  } 


T-t 

-T 


(IX.  13) 


+  U(t>2T)^ 

-T  -T-t  T-t 

where  we  have  separated  the  cases    where  the  time  of  interaction  is 

longer  or  shorter  than  the  time  2T  required  to  cross  the  collision 

sphere. 

In  a  similar  manner  we  evaluate  G^  distinguishing  between  the 

following  cases: 

Case  I:     -T<t    ;t    +  t  <-t 
o  o 


or   T<t;t  +t<T 
o  o 


g  (t,  t^,  p,v)=  g^(t,  t^,  p,v)  (IX.  14a) 


Case  II:   t  <  t   <  T;  T  <  t    +  t 
o  o 


g  (t,  t  ,p,v)=  g  (t,  t  ,p,v)  (IX,  14b) 

O  w  o 
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Case  III:      t   <-T;-T<t    +  t<-T 
o  o 


g  (t,  t  ,  p,  v)  =  g    (t,  t  ,  p,  v) 

O  JO 


(IX.  14c) 


Collisions  which  enter  the  strong  collision  sphere  are  neglected  because 
o£  the  strong  oscillations.     This  yields 


Gjt,  p,v)=  U(T-T>t) 

D 


r-T-t     rT-t       fT   '  f-T 

U.dt    +|f    dt    +1    f   dt    +1  $  dt 
flollol  2ol3o 

J-rp  JX-t  •'T  t 


T-t 

,    .  {IX. 15) 

T-t 


+  U  (t>T-T) 

where  again  interaction  times  longer  or  shorter  than  (T-t)  have  been 


separated.    In  the  expressions  for  G    and  G    we  realize  after  a  change 

a  b 

of  variables  that  the  corresponding  integrals  over  $      and  $  are 
identical.    From  the  Eqs.  (IX.  11a)  and  (IX.  12)  it  is  also  clear  that 
^  ^    is  a  symmetric  functi( 
integral  one  finally  obtains 

i  -G  (t,p,  v)  =  U(2T>t) 


is  a  symmetric  function  in  z  =  t^  +        .     Performing  the  1^- 


^T-t         \  (IX.  16a) 


[  r$,dt  +r  $^  dt  \ 

>l   1    o    I      2  o( 

^^;*'-t/2     •'T-t  ^ 

+  U  (t>2T)|J l^^^o"^  ^"T"  '^^  '  ^  4I 


and 


h'  G^(t,  p,  v)  = 


=    U(T-T>t)|  r$,dt    +A-dt  \ 

If 


+  U(t>T-T)<;  l^^dt^  S     .  (IX.  16b) 
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We  now  introduce  the  following  dimensionless  variables 


D 


_o 

"o  D 


J- 

1^ 


kT 


4un  e 
e 


s  =      t  with  uj    =  */2' 

D  P  V 


UU  = 
P 


Srrn  e 
e 

m 


y  =  t_ 


,/T 


u  =   —       with  V 

V  av 
av 


1  m 


(IX.  17) 


and  the  following  abbreviations 


R  = 


u  •  s 


P  =   -r:  = 


T 


X     -  X 


(IX.  18) 


With  these  definitions  the  preceding  relations  can  be  rewritten  as 


gj(s,y,x,u)= 


2  2 
X    +  (1-x  )y  (y+R) 


2  2,  2 

X  +  (1-x  )y 


2 

x  +  (1-x 


)(y+R)^ 


g2(s,  y,  X,  u) 


g^(s,  y,  X,  u) 


xu 


1- 


X    +  y  (1-x  ) 


V-2  r-  2" 
X    +  y    (1  -X  ) 


(IX.  19) 


2        r— 2^ 
=    —  •  -t/  1  -X 
xu      f  40 


and 


^(s,  y,  X,  u)  =  cos  I  C-  g^{s,  y,  x,  u)  |  -1  (JX.  20) 


where 


C  =  C  •  C  (IX.  21a) 


and 


C    =     4"  -n'q')  (IX.  21b) 

1         Z        c  c 


C   =  ^ 


^  =  0.  03043t/-^i^  .  J^t^ 
2kT  Y        18  T 


2        m- D- V  2kT  _     ,  ^. 

av  '10 


Similarly  we  have  for  the  integrals  over  t^ 


G^(s,  X,  u)  =  U(2>R)|  r$  ^dy  +[  ^^dy  | 

'i^        ^1-R  (IX.  22a) 


+  U(R>2)|J'^$^dy  +  (  ^  -  1) 


and 


G^(s,  X,  u)  =  U(1>R+P) 


1/  ^^"^^  "^^1 

{a  ay} 


+  U(R+P  >  1)<  ^  $^  dy>  (IX.  22b) 


which  leads  to  the  thermal  average 

2rl 


F(s,  n^,  T)  =  2TTn^D  /   du  =,u  •  e "  ^  /     dxx-^l-x    2-G(s,x.  u)      (IX.  23) 

•'o  ^ 
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with 


,     G(s,  X,  u)  =  U(x>x  )  .  G  (s,  X,  u)  +  U  (x  >x)-  G,  (s,  x,  u).  (IX.  24) 

o         a  o  b 

These  integrals  have  been  evaluated  numerically  using  the  program 

(FORTRAN  IV)  discussed  in  Appendix  A.    This  program  calculates  the 

thermal  average  F  as  a  function  of  the  normalized  time    s   for  the 

parameters  (nq  -n'q'),  n    and  T.    The  upper  cutoff  is  given  in  units  of 
c         c  e 

the  Debye  length  and  the  lower  cutoff  in  units  of  the  strong  collision 

cutoff  p     of  Eq.  (IX.  5). 

o 

Before  we  discuss  the  methods  for  obtaining  the  Fourier  transform 
of  F(t)  and  the  actual  intensity  profile,  it  is  useful  to  derive  the  small 
and  large  time  limit  of  F(t).    The  small  time  limit  is  determined  by  the 
integrals  over  f  ^  and  gives  the  asymptote  of  the  thermal  average  for 
the  static  wing.    The  large  time  limit  depends  only  on  the  $  ^  integrals 
and  yields  the  thermal  average  as  required  by  the  iinpact  theory. 
In  the  small  time  limit  $  ,    reduces  to  the  form 

-1 

(IX.  25) 
where 

^  +  v^t^      .  (IX.  26) 

o 

This  expression  depends  only  on  the  instantaneous  distance   r  as 
expected  in  the  static  limit  and  the  thermal  average  is  therefore 
obtained  immediately  by  the  integral  over  r 

f^t),  ^=4nn    pr^$  (t.  r)       dr     .  (IX.  27) 

t-»0  el  1  t-»0 

•'o 


=  cos 


—  (nq  -n'q')  j 

2         c         c    m  2 
r 
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3  h  t 

In  the  small  time  limit  where   — (nq  -n'q')  —  '  — -       0     we  can 

2       c         c     m  Z 

then  perform  the  integral  with  the  result  c 

2    r  1— 

n  =  "  pMrn  -n'q'   )  —  t    2  (IX.  28) 

t-*(J  3     el  c  c     m  J 

For  the  limit  of  large  times  of  interaction  we  have  to  solve  the 
integral 

p(t)       =  2^n  I  dv  vf(v)|   dp  p  t  .  .  (IX.  29) 


o 


For  simplicity  we  set   p^  equal  to  zero  (for  p        0  see  Appendix  B). 
After  a  change  of  variables  and  a  partial  integration  the  integral  can 
be  rewritten  as 


r   4  2  -u^  T  "°  j 

I  du  ==  u  e        CI   —  — 


F(t)       =  -  2nn  tD^v     /  du        u"e   "    Cl   —  ^  '  dz 

t-^  e  avf 


(IX.  30) 


The  z  -integral  is  known  as  Raabe's  integral  (see  p.  144  of  Bateman,  1953) 
and  can  be  expressed  in  terms  of  exponential  integrals.  Furthermore, 
from  Eq.  (IX.  21)  we  realize  that  for  most  practical  situations  C  <<1. 
Keeping  only  the  leading  term  in  C  we  have 


(IX.  31) 
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where 

2        2nn^      /  3(nq^-n'q'he  \2 
4C    =     -7^^  Tt~"^  J    and  B=  0.27  (IX.  32) 

The  large  time  limit  of  the  thermal  average  in  Eq.  (IX.  31)  is  required 
for  the  calculation  of  the  line  center  and  all  modern  impact  theories  give 

the  same  result  except  for  the  additive  constant   B  whose  value  de- 
pends on  the  particular  cutoff  procedure  applied.    Appendix  B  gives 
a  summary  of  the  different  constants  obtained  in  the  literature  which 
vary  considerably.     To  what  extent  this  uncertainty  shows  up  in  the 

final  line  profile  depends  on  the  value  of  the  constant  C.     The  influence 

2 

will  be  small  if  -  In  (4C  )  is  considerably  larger  than  the  uncertainty 
in  the  additive  constant  B.    Furthermore,  the  large  time  limit  of  the 
thermal  average  affects  primarily  the  center  of  the  line  profile  and 
its  contribution  vanishes  when  moving  into  the  line  wings. 

Finally  we  show  numerical  results  for  F(t)  as  obtained  by  means 
of  the  program  in  Appendix  A.  Most  of  the  calculations  shown  in  this  paper 
have  been  performed  for  the  following  electron  density  and  temperature 
parameters. 


case 

n^[cm-^] 

^e  [^] 

experiment 

A 

8.  4-  10^^ 

12200 

Boldt  and  W.  S.  Cooper,  1964  (cascade  arc) 

B 

3.  6.  10^"^ 

20400 

R. 

C.  Elton  and  H.  R.  Griem,  1964(T-shock  tube) 

C 

1.3.10^^ 

1850 

C. 

R.  Vidal,  1964,  1965  (RF-discharge) 
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These  parameters  correspond  to  experiments  which,  as  stated  already 

in  the  introduction,  have  revealed  the  largest  discrepancies  between 

experiment  and  the  modified  impact  theory.    We  will  concentrate  our 

attention  on  the  high  density  case  A  and  the  low  density  case  C,  since 

case  B  is  regarded  as  being  less  accurate  because  of  lacking  absolute 

intensity  calibrations. 

Figures  2  and  3  show  the  normalized  thermal  average  F/  F  as 

a  function  of  the  dimensionless  variable  s  =  of  .  t  for  the  cases  A  and  C. 

P 

Figure  3  shows  the  results  for  three  different  Stark  components 

specified  by  the  quantum  numbers  n    =  nq  -n'q'  •  F     is  the  small 

k  ceo 

time  limit  according  to  Eq.  {IX.  28)  whose  Fourier  transform  leads  to 

the  static  wing.    The  dashed  lines  are  obtained  with  a  lower  cutoff 
2 

n     .    =  0    ='^+na.    It  can  be  seen  that   for  case  C  the  dashed  curves 
mm       o  o 

get  closer  to  the  static  limit  F^  than  for  case  A.    In  order  to  obtain 

the  thermal  average  F   for  the  limit  n     .      0  the  numerical  calculations 

mm 

were  finally  performed  with  typically   n     .    =^  0.  01  n     so  that  F     ,  and 
^  ^  y    ^mm  ^o  calc 

F    differed  less  than  about  0.  1%  over  at  least  one  order  of  magnitude 
o 

in  s.    Far  smaller  values  of  s,  where  F     ,    and  F    start  to  differ 

calc  o 

again,  F         is  then  replaced  by  F  .    In  this  maiuier  we  obtain  the  solid 
calc  o 

curves  in  Eqs.  2  and  3  which  are  used  in  the  following. 

It  should  be  noted  that  these  curves  are  calculated  on  the  basis 

of  the  dipole  approximation.    It  is  clear  that  for  impact  parameters 
^  2 

p  <  n  a      higher  multipole  terms  have  to  be  considered.    Since  the 

values    s    of  interest  are  approximately  given  by  s  <  of  /ao),  one  expects 

P 

higher  multipole  terms  to  be  less  important  the  closer  F  gets 

_  calc 
— 2 

to  F    for  p        =  n  a   .    This  is    consistent  with  the  experimental  fact 
o  min  o 

that  in  case  A  an  asymmetry  of  the  line  has  been  observed  which  cannot 
be  explained  within  the  dipole  approximation,  while  in  case  C  no 
asymmetry  has  been  observed. 
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For  large  s    Figs.  2  and  3  show  the  transition  to  F    as  given 

00 

in  Eq.  (IX.  31),  which  forms  the  basis  for  the  familiar  impact  theories, 
X.     THE  FOURIER  TRANSFORM  OF  THE  THERMAL  AVERAGE 

Having  calculated  the  thermal  average  F(t)  we  now  focus  our 
attention  on  the  evaluation  of  its  Fourier  transform 


i(Aa)j^)  =— /     exp  (iAojj^s)  F{s)  ds  (X.l) 


as  required  by  Eq.  (V.  17)  (see  also  Eq.  (IV.  I6))where  the  dimension- 
less  variable 

A'D^  =  (AUD-AU).-  e)/(ju  (X.  2) 

Rip 

is  the  frequency  separation  from  a  particular  Stark  component 
(cf.  Eq.  (V.  14))  for  an  ion  field  strength  6  in  units  of  the  plasma 
frequency  uu   .  . 

The  thermal  average  F(s)  does  not  immediately  allow  a  straight- 
forward  Fourier  transformation  because  for  large  s  F(s)  is  proportio- 
nal to  s  according  to  Eq.  (IX.  31),  hence  i{/\UJ-^)  diverges.    This  diver- 

R 

gence  is  due  to  the  fact  that  we  neglected  the  finite  lifetimes  of  the  un- 
perturbed states  involved  which  naturally  terminate  the  maximum  time 
of  interaction  s.     This  may  be  taken  care  of  by  introducing  a  convergence 
factor  exp  (-gs)  which  can  be  obtained  by  replacing  the  delta  function 
in  the  power  spectrum  of  Eq.  (3)  in  paper  I  by  a  narrow  Lorentzian 
line  with  a  natural  width  e   (E.  Smith  and  Hooper,  1967).    In  the  final 
line  profile,  however,  natural  line  broadening  is  always  negligible  with 
respect  to  Stark  broadening  which  allows  us  to  set   e  to  zero  without 
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affecting  the  shape  of  the  profile.    For  this  reason  we  will  evaluate 


On/ 


i(AuUj^)  =^lim  ^  /      e  "^^  e  F(s)  ds  (X.  3) 


F(s)  is  known  numerically  and  there  are  many  ways  to  perform  the 
Fourier  transform.    In  order  to  find  the  most  convenient  method  we 
notice  that  according  to  Eqs.  (IX.  Z8)  and  (IX.  31),  'F(s)  has  the  following 
asymptotes 


_  3/2 
for   s-0:     F  (s)  =  p  s  ' 
o  1 


and  for         s-»oo:     F  (s)  =  p  s  , 

00  2 


{X.4) 


where 


and 


2  3  3/2 

p    =  -  — n  D    (2nC)  ' (X.  5) 
i         o  e 


=  -  4^n^  [b  -  In  (4C^]  , 

The  transition  from  F    to  F    is  very  smooth  because  the  power  in  s 

O  00 

changes  only  by  l/2  over  the  entire  range.    It  has  been  found  that  F(s) 
may  be  approximated  by  a  function  G(s)  whose  Fourier  transform  can 
be  given  analytically  and  whose  parameters  may  be  determined  by  a 
least  square  fit.    The  function  G(s)  can  be  given  in  terms  of  the  series 

G(s)=^Gj^(s)   ,  (X.  6) 

k 

where  the  number  of  terms  in  the  series  depends  on  the  required 
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accuracy.    As  a  first  approximation  Eq.  (X.  4)  suggests 

2 

a^s 


G  (s)=    .  ^  ,  (X.  7) 


2b^s 


with 


aj_  =  P2  (X.  8) 


and 

bj  =  i(p^/Pj)2 


Gj^(s)  has  the  small  and  large  s  behavior  of  F(s).  It  then  turns  out 
that 

for    s-0      F^(s)-G^(s)  =  p^s^^^ 


and  for    s^oo       F(s)  -  Gj(s)  =  p^   ,  (X.  9) 


where  p^  and  p^  now  have  to  be  determined  numerically.  Consequently 
we  take  G  (s)  to  be 

5 

3.  S 

°2'«'=   ,  2  '       ,5/2      •  '^-l"* 
(s  ^-Zb^s)  ' 

It  then  becomes  apparent  that  G  (s)  is  given  by 

3k-l 

=        ,  2^^!;     ,21c- 3/2  '^•"> 

(s  +2b,  s)  ' 
k 
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4k-  3 


such  that  one  obtains 


for   s-»0:      G{s)  =  y    p^^  ^  s  ^ 


00 

2-k 


—  K 

k=l 

In  this  manner  the  Fourier  transform  of  any  G  (s)  can  be  expressed 

in  terms  of  modified  Bessel  functions    K    and  K, .    For  all  situations 

o  1 

calculated  it  was  found  that  G  (s)  and  G  (s)  were  sufficient  to  keep 
the  deviation   F(s)-G{s)  smaller  than  1%  for  all  values  of  s.    In  some 
situations  a  fit  better  than  2%  was  obtained  with  Gj^(s)  alone.    As  a 
further  advantage  it  should  be  noted  that  this  method  tends  to  suppress 
"noise"  introduced  by  the  numerical  evaluation  of  F(s). 

In  the  following  we  evaluate  the  Fourier  transform  i(k,  AuJj^)  of 
any  G  (s)  as  defined  by 


1     r  -es  ^^^R^ 
i(k,A«;j^)=  Hm  —   /    e        e  G^(s)  ds  (x.  14) 

•'O 

Their  sum  will  then  give  us  the  desired  Fourier  transform  i  (Aa)-o)'  In 
particular  we  are  interested  in  i(k=l,Au)„)  and  i{k  =  2,  Auu^^).    We  have 


2 

xAuu  ^  s  as 
-es     "  R  1 


ws  +2b  s 
o  If  1 


J  iAti)  „  s 

1     r  « 
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Introducing 


(X.  16) 


one  finally  obtains 

2 


=  a^-b^    (cos  Z^-  i  sin  Z^)  (X.  17) 


where  H^^^and  H^}^  are  Hankel  functions  and  J  ,  J  ,    Y    and  Y  Bessel 
o  1  o      1        o  1 

functions.    These  functions  like  all  the  other  functions  used  in  this 
report  are  consistent  with  the  definitions  as  given,  for  example,  by 
the  NBS  Handbook  of  Mathematical  Functions  (Abramowitz,  1969). 
For  large  arguments  is  also  useful  to  have  the  asymptotic 

expansion 

2 


•  /        X         3      ^1*1  1    (i   L    .i.   5         .,   3.  5-  7  ., 


9-  25-  7  9-  25-  49*  11 

-d-i)  +.   —   (1+i)  +  -  - 


2-8  Zl  8^Z^  )  (X.18) 

.  LJT,  ,  L      _L  (1^)    .  WM  W 

8t  ,  Pi     R       I  4   ^pj     A.^         32  \p^; 


2 

2/  AUJj^ 
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Using  Eq.  (X.  5)  for  p    the  latter  relation  gives  us  exactly  the  Holtsmark 

-5/2  ^ 

AX        wing  for  all  Stark  components 

yAuOj^)  =  nn^D^C^/^-  Aa)"j^^^   .  (X.  19) 


In  a  similar  way  one  derives 

1  -es    ^^"^R^  ^2^" 

^(AuUt,)  =  i(k=2,  Aa)„)  =    lim — 4  e,        e  — 5  E~7^  ds 

^      ^  R  rrjf^  (s^2b2s)^/2 


(X.  20) 


With 


^2  "  V^R  (X.21) 


one  finally  obtains 

^2^2     '^h  \  n)  2 


} 


(1 )  2 
+  H^^  '{Z^)  {IGZ^  +  1282^-3) 

a^b^  (X.22) 
=  — ^         (cos  Z^-i  sin  Z^) 


if- 


36Z  J  (Z  )  +  J.(Z  )  (l6Z^-3)  -  Y  (Z  )  (16Z^-15)  -  28Z  Y  (Z  ) 
11         2o2         12  2  o2  2  2-12 


+  i  I  JJZ^)  (16Z^-15)  +  28Z2Jj(Z2)  -36Z2  Y^(Z2)+Y^(Z2)(l6Z^-3) 
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The  asymptotic  expansion  for  large         is  given  by 


128Z^ 


) 


(X.  23) 


If  one  requires  an  even  better  fit  of  G(s)  to  F(s)  the  general  transform 

i(k,  A«J    )        defined  in  Eq.  (X  .  14  )  is  given  by 
R 


4^'02k-2)!  1  .  -  .k+1  d 


(4k-4)! 


a,   lirrx  —  ( - 1 ) 

TT 


k+1    ^2k-2    /  b^(e-iAti3j^) 

e  ^    m  (e-iAt«p) 


k+1  ,2k- 2, 
d      e    d  b 


j^b^(e-iAt«j^)| 


"o  I  "k^"  "  R 
(X.  24) 

Finally  we  want  to  show  that  this  technique  always  gives  the  static  wing 
according  to  Eq.  (X.  19)  for  large  Atu.    For  this  purpose  one  has  to 
perform  the  Fourier  transform  of  the  small  time  limit  of  G{s)  as  given 
in  Eq.  (X.  13). 


-^""-R'  -  Air-^oo  ^^^-R'  =  2  P2k-1 


on  U  A,.,  »L       n  ' 

R 


k=l 


e       e  s  "^ds 


(Zk+1)! 

2k+l  P2k-1 
2  k! 


exp 


i-^{4+k)l 


2^2'  "'J  ^^R 


-(k+-) 


(X,  25) 


5  P3 


2    Auu„  4  2 

R  AoJj^ 


5- 7- 9  7 


8 


(1-i)  -  + 


R 
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One  recognizes  that  the  first  two  terms  are  identical  with  the  first 
terms  in  the  Eqs.  (X.  18)  and  (X.  23).  Hence,  we  always  obtain  the 
static  wing  for  large  AuUj^. 

Another  important  property  of  i(A(JUT-,)  is  that  for  small  its 
leading  terms  in  the  expansion  are 

^    1    .  <^N-^2' 


M.v%^i--o«-K>=-^;;;r-  ....  ^  ^^^^^ 


j  In  this  manner  it  smoothly  goes  over  to  the  Lorentz  profile  of  the  un- 

j  modified  impact  theory. 

Before  discussing  the  numerical  results  of  i(AUU^)  we  first  list 

R 

■I  the  constants  a,  and  b,  for  the  cases  A,  B  and  C  as  specified  at  the  end 
k  k 

I  of  Sec.  IX.    a    and  b    are  determined  from  Eq.  (X.  8),  where  p,  is 
11  1 

1  given  by  Eq.  (X.  5)  and  p^  is  taken  from  the  large  time  limit  of  the 

computed  F(s).     p  as  calculated  by  the  program  of  Appendix  A 

2  comp. 

may  differ  slightly  from  p^  as  defined  in  Eq.  (X.  5),  if  C  is  not  very 

much  smaller  than  unity  because  p^  is  based  on  Eq.  (A.  18), 

2  comp. 

which  is  true  for  any  value  of  C  and  goes  over  to  Eq.  (X.  5)  for  small 

C.    a    and  b    are  determined  numerically  by  a  least  squares  fit. 
If  — 

'    The  maximum  deviations  from  F(s)  obtained  with  G^(s)  alone  and  with 
k  Gj(s)  +  G^ls)  are  listed  too. 
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In  presenting  the  numerical  results  of  i(AtD^)  we  concentrate  on  the  real 

part  which  turns  out  to  be  the  most  important  part.    We  have  chosen 

two  different  normalizations.    In  Figs.  4  and  5,  i(Auu    )  is  normalized 

R 

with  respect  to  the  large  frequency  limit  i  {/Sw^)  to  show  the  useful 

CO  K 

range  of  the  static  theory.     The  short  vertical  lines  mark  the  position 
of  the  Weiskopf  frequency 

Ati)    =  V    V  f4-<^^  -^'q')  ~l  =  ^  /C  (X.  27) 

c       av  '   \  2       c         c    m/       p'  ^ 

for  a  particular  component  (nq^-n'q^)  which  according  to  classical 

arguments  determines  roughly  the  range  of  validity  for  the  static 

theory  (see  p.  321  of  Unsold,   1955  and  paper  II).    It  should  be  pointed 

out  that  Auu     is  usually  defined  in  terms  of  an  average  Stark  splitting.  In 
c 

both  cases  A  and  C  /\uu    describes  the  range  of  the  static  theory  very 

well.    If  one  allows  for  a  deviation  of  about  10%  at  the  most  from  the 

static  asymptote,  Aoj^  may  be  lowered  effectively  by  more  than  an  order 

of  magnitude.    A  more  detailed  discussion  is  given  later  with  the  final 

line  profile  calculations. 

The  other  normalization  with  respect  to  the  small  frequency 

limit   i  (Auu^)  is  shown  in  Figs.   6  and  7  for  cases  A  and  C  again.  These 
o      R  °  ° 

plots  demonstrate  the  useful  range  of  the  unmodified  impact  theory, 

which  is  based  on  i  (AuUt^)  and  is  expected  to  break  down  around  the 

o  R 

plasma  frequency,   as  can  be  seen  in  Figs.   6  and  7.    In  order  to  extend 
the  range  of  validity,  the  modified  impact  theory  makes  an  impact  para- 
meter cutoff  at  v/auu  (  the  Lewis  cutoff)  whenever  this  is  smaller  than 
the  Debye  length  D;  this  cutoff  accounts  for  the  finite  time  of  interaction 
to  second  order.    More  details  are  given  in  Appendix  B.    The  corres- 
ponding function   i,       .    (Auu^)  has  been  included  in  Figs.   6  and  7.  Since 

Lewis  R 

the  usual  derivation  of  i^       .  (AO)-.)  is  based  on  the  limit  of  very  small  C, 

Lewis  R 

one  expects  the  best  agreement  between  the  Lewis  result  and  our  result, 
which  considers  the  finite  time  of  interaction  to  all  orders,  for  the 
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situation  with  the  smallest  C.     That  this  is  in  fact  true  can  be  seen 

from  the  low  density  case  with  nq^-n'q^  =  3.     This  component  is  plotted 

again  in  Fig.   8,   in  order  to  demonstrate  the  importance  of  G  (s)  for 

those  cases  where  the  deviation  of  G^(s)  from  F(s)  is  large  (Table  I 

gives  a  maximum  deviation  of  13%). 

Figures  6  and  7   also  contain   the  static  limit  i  (Auu„)  (dashed 

00  R 

lines)  and  the  Weiskopf  frequency  /\w   .    It  gives  an  idea  how  close  the 

Lewis  results  get  to  the  static  limit.     One  notices  that  with  increasing 

values  of  C  the  deviation  of  i^       .  (Auu„)  from  the  static  limit  becomes 

Lewis  R 

larger.    In  his  line  wing  calculations  (Griem,  196Z,  1967a)Griem  adjusts 

his  "strong  collision  term"  E  ^  ,  in  such  a  manner  that  the  Lewis  result 

BP 

is  identical  with  the  static  limit  at  the  Weiskopf  frequency.    In  the  Figs. 

6  and  7  this  means  that  the  straight  line  representing  i^       .    (auu^)  is 

Lewis  R 

shifted  to  the  right  until  it  cuts  Auu   .    We  use  here  Au)    as  defined  in 

c  c 

Eq.  (X.  27)  for  every  individual  component  instead  of  the  average  value 

A^^  -  kT/(hn  )  used  by  Griem.    Since  the  Lewis  line  would  then  lie 

appreciably  above  the  curve    i(Auo    )  one  realizes  that  this  procedure 

R 

definitely  overestimates  the  electron  broadening  as  already  observed 

experimentally  (Vidal,   1965;  see  also  Pfennig,   Trefftz  and  Vidal,  1966). 

A  better    method  would  have  been  to  adiust  E  ^  ,  such  that  i^       .  (Auu„) 

Bp  Lewis  R 

forms  a  tangent  of  the  static  limit.    However,   it  is  clear  that  any 
adjustment  of  -^g^/  effectively  changes  the  range  of  the  unmodified 
impact  theory  and  also  defeats  the  purpose  of  the  Lewis  cutoff,  namely 
to  correct  the  completed  collision  assumption  to  second  order. 

Finally  it  ought  to  be  emphasized  again  that  except  for  the  time 

ordering  the  Fourier  transform  of  the  thermal  average  i(A«3T3)  ^®  pre- 

R 

sented  here  takes  into  account  the  finite  time  of  interaction  to  all  orders. 

Hence,  for  small  Aw_  it  goes  over  to  the  impact  theory  limit  and  for 

R 

large  Ati^-o  it  gives  the  static  limit  without  requiring  a  Lewis  cutoff. 
R 
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XI.     THE  ONE-ELECTRON  LIMIT  FOR  HYDiROGEN  AND  THE 
ASYMPTOTIC  WING  EXPANSION 
Having  obtained  the  Fourier  transform  of  the  thermal  average 

i(AuJTD)  now  prepared  to  calculate  the  actual  line  intensity  by 

R 

evaluating  1(a), according  to  Eq.  (IV.  15)  and  averaging  it  over  all 
ion  fields  according  to  Eq.  (II.  1).    As  explained  in  Sec.  IV.  this  problem 
is  greatly  simplified  in  the  one  electron  limit  where  no  matrix  in- 
version is  required  and  the  intensity  I(A(Ju)  is  given  by 

/•oo 

I(A(Ju)  =  UAoj)  +  1    P(e)  I  (Am,  B)  dQ    .  (XLl) 


I(Auu)  is  the  static  ion  contribution  originating  from  the  first  term 
i 

1/Auu^P  in  Eq.  (IV.  21 )  and  I(A(D,  B)  is  given  by 

I(AU0,S)=    —  y^(nq  m  jdln'q   '  m '  )  (  n 'q '  m '  I  dl  nq  m  ^ 
Tj        ^     aa  aa  bb  bb 


(XL  2) 

dt  expfiAtu^tJ  (n'q'm';  nq  m  l'^^^  ^(t)|  n'q'm';  nq  m  > 

bbbb  aaaa 


using  the  definitions  of  Eqs.  (V.  14)  and  (X.  2).     The  density  matrix  p 
is  assumed  to  be  constant  over  the  relevant  initial  states.    With  Eq. 
(IX.  1)  the  last  expression  can  be  rewritten  as 

I(Atu,P)=  V(n;   m  Idln';/m')(n'q'mMn';/m'>(n'/, 'm'ldlrtf,  m  ) 
'  aa  a    a  bb        cb  cb  cb 

(n^  m,  Ina  m,  )  (n';, 'm' I  n'q'm')  (n'q'm' I  n'/  'm')<f  n'^ 'mM  n'q'm'^ 
cbbb  ac         cc  cc         be  bb  bb 

-m  -m, 

( nq,  m,  I  nJ  ,  m,  )  <  nJ^^,  m  Inqm)(nqm  Irtf,  m)(-l)  (2L+1) 
b    b       b    b"^       b    c       c    c         c    c       a  c 

'I'  r    L  \ /;/  ^    L  \ /;  '  i.  L  \  /  '  i  L\ 
aa       \/aa        \/bb  \/bb\ 

li(Aa.j^,6,n,n',  q^,  q;>q^.q^) 

 „      'm'm  M/\-m'm  M/\-m/m,  M/ 

c    c        ^^t^    ^  (XI.  3) 


where  the  dipole  matrix  elements  have  been  transformed  from  parabolic 
to  spherical  states  and  the  summation  over  intermediate  states  I  nq  m  ) 


a  a 


and  In'q'm')  has  been  performed.    We  next  apply  the  Wigner  Eckart 
theorem  (see  Eq.  (  5.  4.  1 )  of  Edmonds,   I960)  to  the  dipole  matrix 
elements  and  replace  the  reduced  matrix  elements  by  the  corresponding 
radial  matrix  elements  (see  Bethe  and  Salp.eter  1957). 

/ 1    1    I \/l  1  I  \ 

(n^.mld  In';/m')  =  ( - 1  )""^  ij(2^+l  )(2^ '  +  1 1  l|  )<n^,  |  r!  n'^'> 

.  Vm|.i  myyO  0  0  y 

(XI.  4) 

Inserting  this  relation  into  Eq.  (XI.  3)  and  using  the  orthogonality 
properties  of  the  3j- symbols  we  have 

I(Auu,  B)  =  "V^  (n'q/m/  I  n' I 'm/)  <nf,  m^  |  na  m  N  (n'^  'mM  n'q'm') 
/  J         bb        cb        cbbb  ac  cc 

(n'q^mM  n 'i;m^)  (n't^m^l  n'q^m;^  (nq^m^l  m^m^Vnt^m  J  nq^m^) 


(nq  m  Inf   m  >    (2;,  +1  )(2|/  +  1  )(2^  +1  )(2r '  +  1 )  j     <  n^.   !  r  I  n '  )  <  n '  I  r  I  nl,  ) 
ccac    i       a  a  c  cl  a  a  c  c 


\p  0  0/  \  0  0  0/  \-m'  m  MV  \-m'  m  M7  \  -m/  m,  M/ \ -m/m  Mi 
^  ^cc/\cc/\bb/^bb 

i(Auu^,  8,  n,  n',  q^,  q'   q  ,  q'  )  (XI.  5) 

R  b     b     c  c 
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;  J£  we  finally  replace  the  unitary  transformation  by  the  corresponding 
3j-symbols  according  to  Eq.  (VIII.  8)  the  result  is 

I(AaJ,e)=  V(2^  +1)(2^'  +  1)(2;   +1)(2;  '  +  1)(2^  +l)(2f'  +  l) 
jL^      a  a  b  b  c  c 


a     a     \/   c     c     \l    a     a      \/  b     b      \/  b     b      \/   c  c 


0  0/ 


0     0     0/  \  0    0     0/\-m'  m  M7\-m'm  ]w\  -m/  m^M/Vrn/  M> 

cc'^      cc/^  bb/^bb 


n-1 

2 

2 

n-1 

2 

m  -q 

c  c 

m  +  q 
c 

n^l 


y~-  ^b 


m  -q      m  +q 

c     c       c     c  ' 

-r    m  ; 

^             2  ci 


/ 


I' 

c 


2  -'"ci 


b  ^b       b   ^b         ,  M     bjb      b  ^b  , 

 m'  A   —   ; —  -m' 

2  b  /  \      2  2  b 


'b/\  2 
/  \ 


/ 


/  n'-l    n'-l  .  \ 

I    ~1         —  -^a    \  '  2 


m'-q'  m'+q' 
c     c      c  c 


m 


n'-l  ,  \ 


m'-q'  m'+q' 
c     c  c  c 


m 

2  c 


/ 


(np    I  rl  n';  '>(n'^  '  I  rl  nt  ) 
a  a  c  c 


i(^'^R'   ^'   ^'  ^b'   %'   <  ^ 


(XL  6) 
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The  preceding  relations  hold  for  the  general  case  of  upper  and 
lower  state  interactions.     They  simplify  considerably  if  there  is  no 
lower  state  interaction  (e.  g.  ,   Lyman  lines).     Then  one  obtains 

^     I     ^      ^     n-1  n-1 

MAuo,R)=  l(nllrllO)l^y^ 


n-1 
2 


m  -q   m  +q 
c     c     c  c 


 m 

L  < 


i  (Ao)^,  B,  n,  q  q  ) 
u      -R-  be 


(XI.  7) 


with 


i  (AuUp,  B,  n,  q     q  )  =  \Km^,  B,  n,  n'  =  1,   q     q '  =  0,   q  ,  q'  =  0). 
uK  DC  K  bb  cc 


(XI.  8) 


Equation  (XI.  7)  may  be  further  simplified  by  evaluating  the  3j- symbols 

and  summing  over    m,  and  m    with  the  result. 
^  be 


I  (Am, 8) = 
u 


I(nl|rll0)| 
4(n^-l)^n^ 


n-1 

■E 


q  +n 

2  ^22 
n  +(-1)         (n  -2q  ) 

b 


n-1 


q  +n 

2  *^        2  2 

n+(-l)  (n''-2q^) 

c 


q^=-(n-l)L 


i  (Auu„,  B,  n,  q  q  ) 
u       XV  be 


(XI.  9) 
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These  simplified  relations  may  also  be  used  for  the  higher  series 
members  of  the  other  series,  whose  transitions  do  not  end  on  the 
ground  state  if  lower  state  interactions  contribute  only  a  negligible 
amount  of  broadening  to  the  final  line  profile. 

The  foregoing  relations  for  the  one  electron  limit  essentially 
represent  the  asymptotic  expression  for  the  intensity  in  the  line  wings. 
If  one  is  interested  in  frequency  perturbations  Auu  which  are  significantly 
larger  than  the  average  ion  field  splitting  Eq.  (XI.  1)  can  be  simplified 
by  replacing  the  ion  field  average  of  the  electron  contribution  by  the 
electron  contribution  for  the  average  ion  field  6 

av 

KAud)  =  I.  (AO))  +  I  (Am,  B^^)  (XI.  10) 

with 

OS 

0    =  /  BW(e)  de  .  (XI.  11) 

av  J 

o 

K  A'i)  is  very  much  larger  than  the  average  ion  field  splitting,  then 

according  to  Eq.  (X.  2)  Auu^,     Am/uu      and  I(A(U,  may  be  replaced  by 

°  rl  p  av 

I(A(JU,  B=  0). 

KAo))  =  I.(A(u)  +  I  (Aud,  P=  0)  (XL  12) 

In  the  limit  B      0  the  Eqs.  (XI.  5)  to  (XI.  95  simplify  drastically  because 
i(A'i)j^)  depends  no  longer  on  the  quantum  numbers  q^  and  q^  which 
specify  the  Stark  components  shifted  by  the  quasistatic  ion  fields.  This 
allows  us  to  sum  in  Eq.  (XL  7)  over  q^  and  m^  which  gives  us  for  the 
case  of  no  lower  state  interaction 
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I  im,  R=o)  = 

u 


1  <nll  rl  10) 


qm 


/  n-1      n-1  A 


m-q  m+q 


i  (auu,  B=0,  n,  q) 
u 


i 


{XI.  13) 


(n-1) 


nil  rl  10)        V^/  2  ,  ,,q+n,  2  ^  2,\  .  , 
 r  y       n  +(-1)^      (n  -2q  )     i  (Am,B=0,  n,  q) 

2(n^-l)n        ^\  /  ^ 


2(n  -1) 

q=-(n-l) 

For  the  general  case  of  upper  and  lower  state  interaction  we  can  sum  in 
Eq.  (XI.  6)  over  q^,  q^,  m^,         and  M   and  after  applying  Eq.  (XL  4) 
we  finally  sum  over  the  intermediate  spherical  states  to  obtain 


I(A(JU 


,3=0)  =  ^^  |(nqmldln'q'm')l^  i(Aa),B=0,  n,  n',  q,  q')  , 

(XI.  1 


4) 


m,  m 

How  far  into  the  line  center  the  simplified  relations  (XI.  10)  and  (XI.  12) 
may  be  used,  depends  on  the  required  accuracy.    Numerical  results, 
which  compare  the  asymptotic  wing  expansions  with  the  more  rigorous 
unified  theory  calculations  describing  the  entire  line  profile,  are  given 
at  the  end  of  the  next  section. 

XIL     THE  UNIFIED  THEORY  FOR  HYDROGEN 
In  those  cases,  where  the  entire  line  profile  including  the  line 
center  is  required,  the  line  intensities  have  to  be  calculated  on  the  basis 
of  the  unified  theory.    It  has  to  be  pointed  out  that  in  principal  even  in 
calculating  the  distant  line  wings  the  unified  theory  has  to  be  used  when- 
ever AdJj^   in  Eqs.  (X.  1)  or  (V.  17)  is  no  longer  large  compared  to  unity. 
This  will  happen  in  the  final  integration  over  ion  fields  whenever  8  is 
close  to 


B    =  Aoj/Aoj.  (n,   q      n,'  q') 


(XII.  1) 
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However,  it  was  shown  in  the  last  section,  that  for  large  Auu 
one  may  use  one  of  the  asymptotic  expansions  in  Eq.  (XI.  10)  and  (XI.  12). 

In  the  unified  theory  we  have  to  evaluate  the  following  expression 

I(AUJ,e)  =  —   \^(nq  m  I  d|  n'q'm')  (n'q'm' I  dl  nq,  m,  > 
TT      /  J     aa  aa  bb  bb 


(XII.  2) 


(n'q'm';  nq^m^l    AO)     -£(AllJ     )  In'q'm'; 
bb       bb    lop  op  J  aa 


nq  m  ). 

a  a 


The  matrix  elements  of  Auo      as  defined  in  Eq.  (IV.  14)  are  diagonal  in 

op 

parabolic  states  and  are  given  by 

(n'q'm';  nq  m  I  AcD      I  n  '  q '  m ' ;  nq  m  )  =  Ac«  -  AO)  .(n,  q  ,  n ' ,  q^)  R 
^     ^a    a       a    a  "   op        a    a       a    a  i  a 

(XII.  3) 

where  AOJ .  is  defined  in  Eq.  (V.  14).    The  matrix  elements  of  i!(Aa)      )  are 


given  by 


^n'q'm';  nq  m  I  ^(Aa)     )  I  nq'm';  nq  m  )  = 
b    b       b    b  op  a    a       a  a 


-iTT.(-l) 


•m  -m, 
a  b 


-.2 

j^AaJ-Aa).(n,q^,n',q^)  Rj   ^(2L  +  1 


(n'q'm'In'1'm'Xn't'm'ln'q;  m  '  )  (  n '  q 'm  '  I  n '^.^m^X  n 't^m  M  n 'q^m 
aa        aa  ac         cc  ^ 


(nq  m  1  nt   m  '  )  ^ nq  m  I  n^  m  )  (n;   m  !  nq^  m^) 

aaaaaccc         ccu«-  ^ 
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I'     I  L 
a  a 


\  -m'  m  M' 


I'     I  L 
a  a 


•  m'  m  M 
a  a 


-ro.'  m  M' 


•  m/  m,  M 
b  b 


i  (  Atw  j^,  e,  n,  n,  q^,q 


(XII.  4) 
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Eqs.  (IV.  16),  (IX.  1)  and  (X.  1).    This  relation  simplifies  sig- 
nificantly in  case  of  no  lower  state  interaction  in  which  case  we  need  the 
nnatrix  elements 


(nq.m  li:(Auu     )Inq  m  )  =  -  in     AtJU -Auo  .(n,  q  )  R 
b    b  op        a    a  i        ^  I 


/  -T, — Ti —  ( ^^1."^!, ^  m  Inq  m  >    (nt  m  Inq  m  > 

'^a  b    b       a    a        a    a       a    a  a    c       c    c  J 


^  ,  q  ,  m 
a     c  c 


i  (AuJt3.  B,  n,  q  q  ) 
u      R  be 


(XII.  5) 


Due  to  the  delta  function  the  matrix  of  the  operator  £    is  then  block 
diagonal  in  m,  which  reduces  the  size  of  the  matrices  to  be  inverted 
to  n  X  n  or  (n-1)  x  (n-1)  depending  on  the  quantum  numbers  n  and  m. 
Furthermore,   Eq.  (XII,  2)  simplifies  in  case  of  no  lower  state  interaction. 
After  transforming  the  dipole  matrix  elements  from  parabolic  to  spherical 
states,   applying  the  Wigner-Eckart  theorem  (see  Eq.  (XI.  4))  and  using 
Eq.  (VIII.  8)  one  obtains 


l{m,Q)=  UnllrllO)!^  ^  (-1) 


n+m- 1  - 
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n-1 


m- 


2  2 
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-  m 


n-1 


n-1 
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\  m-q  m+q 
\   b   1 


-  m 


/ 


Im 


(nq.ml    [auu     -Jl(Auu      )|      Inq  m)    .     (XII.  6) 
b        j      op  op  I      '  a 


In  order  to  keep  the  mathematics  simpler  we  concentrate  in  the 

following  on  the  case  of  no  lower  state  interaction,  because  it  covers  the 

experimental  situations  of  case  A  and    B    and  is  also  a  good  approximation 

to  the  higher  Balmer  lines  of  case  C  (see  the  list  of  references  at  the 
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jend  of  Sec.  IX).    Including  lower  state  interactions  means  at  this  stage 
jonly  a  more  extensive  summation  over  3j- symbols  because  the  crucial 
jfunction  i(Ati)j^,  i^),  the  Fourier  transform  of  the  thermal  average,  has 
jalready  been  evaluated  for  the  general  case  of  upper  and  lower  state 
{interactions. 

Using  the  unitary  transformation  of  Eq.  (VIII,  8),  Eq.  (XII.  5) 
imay  be  rewritten  as 


n+m- 1  - 


(nq,ml£(Auu     )Inq  m)  =  -  iTT(-l) 
b  op  a 


|^Atu-/\uu.(n,  q^)  ^+1 ) 


/  n-1         n-1        \    /  n-1 
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m-q  m+q 
a  a 
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-m 
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n-1 


m-q^  m+q^ 
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q  >  0 
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(Amj^,8,n,  ^^^^ 


n-1 


m-q  m+q 
c     c      c  c 


■m 


(XII.  7) 


where  we  have  used  the  fact  that  i  8,  n,  q^  ,  q  )  =  i  (Aa)_,  S,  n,  q,  ,  -q  ) 

uK  bcuK  be 

and  that  i  (AuUt^,  6,  n,  q,  ,  q    =  0)  =  0.    We  also  realize  that 
u      R  be 


<'nq  -mi<C(Auu     )l  nq  -m>  =  (nq,ml£(Auu     )|nq  m) 
b  op        a  b  op  '  a 


(XII.  8) 


A  computer  program  (FORTRAN  IV),  which  evaluates  I(Aud,  B)  according 
to  Eqs.  (XII.  6)  and  (XII.  7)  and  also  performs  the  final  ion  field  average 
according  to  Eq.  (II.  1)  is  presented  in  Appendix  C.    The  ion  microfield 
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distribution  function  employed  is  the  one  given  by  Hooper,  1968a,  1968b, 
which  differs  less  than  about  1%  from  the  values  determined  independ- 
ently by  Pfennig   and  Trefftz,  1966. 


For  the  experimental  parameters  of  case  C,  Figs„  9  and  10 
show  numerical  results  of  I(Auo,  P  =  10)  for  n  =  6  and  n  =  10.     The  fat 
vertical  lines  indicate  the  relative  intensities  and  the  positions  of  the 
Stark  components  for  the  static  field  (ion  field)  3=10  and  it  demon- 
strates the  electron  broadening. 

Figures  11  to  13  show  the  final  line  profiles  I(At«)  after  per- 
forming the  ion  field  average  for  the  experimental  cases  A,  B  and  C 
(see  end  of  Sec.  IX).    As  a  first  result  it  turns  out  that  for  numerical 
accuracies  of  about  2%  it  is  in  all  3  cases  sufficient  to  consider  only 
Gj(t)  meaning  that   i(A(Dj^,  B,  n,  q^,  q^)  may  be  replaced  by 
ij(Aujj^,  6,  n,  q^,  q^)  as  given  in  Eq.  (X.  17).    Although  according  to 
Table  I,  G^(t)  may  differ  from  F(t)  for  some  components  of  case  C 
by  up  to  13%,  it  turns  out  that  after  summing  over  all  Stark  components 
and  averaging  over  ion  fields  this  difference  F(t)  -  Gj(t)  is  apparently 
smeared  out  over  the  entire  line  profile  and  affects  the  final  line 
profile  by  not  more  than  about  2%.     This  is  very  convenient  for 
practical  calculations,  because  it  no  longer  requires  an  extensive 
evaluation  of  the  thermal  average  anymore,  but  for  most  practical 
situations  it  is  sufficient  to  calculate  the  line  intensities  directly  on 
the  basis  of  Gj^(t)  whose  specifying  constants  a^^  and  b^^    are  given 
immediately  by  the  Eqs.  (X.  8)  and  (X.  5). 

This  is  even  more  true  in  view  of  the  fact  that  the  final  line 
profile  is  partially  affected  by  an  uncertainty  in  the  constant  B  as 
defined  in  Eq.  (IX.  31)  or  (B.  19).     As  summarized  in  Table  II  of 
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Appendix  B    its  actual  value  depends  on  the  cutoff  procedure  applied, 


a  problem,  which  has  not  yet  been  solved  satisfactorily.     The  upper 


the  frame  work  of  the  classical  path  theory  (see  also  Chappell,  J. 

Cooper  and  E.  Smith,  1969).     The  lower  cutoff  parameter,  however, 

which  essentially  replaces  the  dynamic  strong  collisions  not  amenable 

in  a  classical  path  theory,   can  only  be  determined  conclusively  from  a 

quantum  mechanical  theory  which  is  also  able  to  handle  strong  collisions 

and  which  does  not  yet  exist.     The  constant  B   adopted  here  is  based  on 

3  2 

a  lower  cutoff  parameter  p     .    =^  +  -r-  n  a  ,  which  specifies  approxi- 

mm  Z  o 

mately  the  region  of  validity  for  the  classical  path  theories  (see  paper 

I).    Numerical  results  based  on  other  values  for  the  constant  B  as 

used  in  the  literature  (see  summary  of  Appendix  B)  are  also  included 

in  Figs.  11  and  12  for  the  cases  A  and  B.     The  largest  value  B  =  1.  27 

is  the  one  adopted  in  the  recent  calculations  of  Kepple  and  Griem, 

3  2 

1968,  while  the  smallest  value  of  B   is  obtained  for  p     .    =     +  -r-n  a 

mm  2  o 

and  choosing  an  upper  cutoff  of  p  ^*  ^06D  as  proposed  by  Chappell, 

J.  Cooper  and  E.  Smith,  1969.    For  case  C   this  variation  of  the  con- 
stant B    does  not  show  up  in  Fig.  13  and  amounts  to  an  intensity  change 
of  at  the  most  about  4%.     These  variations  indicate  the  reliability  of 
the  classical  path  theories  and  demonstrate  that  for  some  cases  the 
error  estimates  given  in  the  literature  are  too  optimistic.     The  effect 
on  the  final  line  profile  due  to  the  uncertainty  of  the  constant   B  will 
be  small  if  either  according  to  Eq.  (IX.  31) 


cutoff  parameter  a  =  p 


/D  (see  Appendix  B)   and  therefore  also  the 
'V  (t')  dt'    can  in  principal  be  decided  within 


-2  In  (2C)  »  1 


(XII.  9) 
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or  if  (Like  for  the  higher  series  members)  the  number  of  Stark  com- 
ponents is  large  which  tends  to  smear  out  the  influence  of  the  constant 
B.    It  should  be  pointed  out  again  that  the  unified  theory  is  intrinsicly 
normalized  independent  of  the  value  of  the  specifying  constants  of  a 
particular  line.    Hence,  any  variation  of  the  constant  B   does  not  affect 
the  normalization  of  the  line  profile. 

In  comparing  the  numerical  results  for  case  C,  with  the 
experiment  it  has  to  be  kept  in  mind  that  we  are  comparing  the  higher 
Balmer  lines  with  calculations  for  the  higher  Lyman  lines,  because 
our  final  line  profile  calculations  have  not  yet  taken  into  account  lower 
state  interactions.     This  means  that  in  a  plot  of  the  intensity  versus  the 
wavenumbers         which  is  essentially  an  energy  scale,  the  line  profiles 
cannot  be  expected  to  coincide  because  of  the  difference  in  the  Stark 
effect.    This  gives  rise  to  different  static  wings  as  explained  in  detail 
by  Vidal,  1965.    Hence,  we  have  to  rescale  the  Lyman  profiles  preserv- 
ing normalization  in  order  to  be  able  to  compare  the  measured  profiles 
of  the  Balmer  lines  with  the  calculated  profiles  of  the  Lyman  lines. 
This  means  that  in  a  plot  of  log  I  versus  log  /\v   we  can  compare  the 
line  shape  of  the  corresponding  lines  directly.     The  agreement  is 
remarkable.    For  the  higher  lines,  n  >  8,  where  Doppler  broadening 
was  shown  to  be  negligible  and  where  lower  state  interactions  no 
longer  affect       the  line  shape  noticeably,  the  agreement  is  better  than 
2%  over  the  entire  measured  line  profile,  which  for  n  =  8  extends  over 

3  orders  of  magnitude  in  intensity.    In  particular,  the  calculations  show 

-5/2 

also  the  surprisingly  large  range  of  the  f\iu         -  wing,  which  extends 
to  l/lO  of  the  maximum  intensity.    This  fact  is  not  explainable  by  a 
purely  static  theory  considering  also  shielding  effects.    For  the  lower 
lines  the  calculated  profiles  have  to  be  folded  into  a  Doppler  profile  in 
order  to  achieve  similarly  good  agreement.    For  the  lower  line  we 
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also  expect  in  the  line  center  some  influence  due  to  lower  state  inter- 
j  actions  on  the  line  shape,  which  is  partially  removed  again  by  Doppler 
broadening. 

-5/2 

A  more  detailed  study  of  the  Aut)         -wings  reveals  some  other 

interesting  facts.    In  Fig.  13,  the  dashed  lines  indicate  the  asymptotic 

5/2  5/2 
A(D         -wings;  except  for  n  =  5  and  n  =  6,  what  appears  to  be  a 

-wing  in  the  measurements  and  calculations  is  not  the  asymptotic 
-5/2 

Holtsmark  /\(U         -wing  in  the  region  of  interest.    If  one  extends  the 

calculations  to  even  larger  frequencies  At)U>  a.11  the  wings  will , eventually 

coincide  with  their  asymptotic  limit.    In  the  paper  of  Vidal,  1965, 

Table  II  gives  a  list  of  the  electron  densities,  which  were  evaluated 

-5/2 

under  the  erroneous  assumption  that  the  measured  Atw         -wing  was 
the  asymptotic  Holtsmark  wing;  it  was  stated  that  for        to  ^-^^ 
electron  densities  coincide  within  i  4%,    A  more  careful  analysis  of 
the  values,  which  have  been  plotted  again  in  Fig.  14  reveals  a  syste- 
matic trend.  For  large  and  very  small  principal  quantum  numbers  the 
electron  density  values  rise  above  the  average  value,  while  the  min- 
imum value  was  obtained  for  n  =  7.    From  Fig.  13,  it  now  becomes 
apparent  that  the  electron  densities  based  on  the  asymptotic 
Holtsmark  wing  will  go  up  for  increasing  n.    For  smaller   n  the 
quantum  number  dependence  of  the  electron  density  is  masked  by 

I  Doppler  broadening  which  raises  the  wings  again  and  explains  the 
increasing  values  of  electron  density  for  small  n.    Another  important 
result  can  be  seen  from  Fig.  13.    For  small  principal  quantum  numbers 
the  line  intensities  are  much  smaller  than  predicted  by  a  quasistatic 

^  theory.    This  was  observed  first  by  Schluter  and  Avila,   1966  and  the 
effective  electron  densities  for  a  quasistatic  theory  as  a  function  of  AX 
show  the  qualitative  behavior  measured  by  them  after  unfolding  the 
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Doppler  broadeningo    This  observation  together  with  the  measurements 

of  Boldt  and  Cooper,  1964^  suggested  the  semiemperical  procedure 

proposed  by  Edmonds,  Schluter  and  Wells,   1967.    A  detailed  quantitative 

comparison  requires  for  the  first  series  members  a  consideration  of 

lower  state  interactions,  which  is  in  process. 

For  the  parameters  of  case  C,  Kepple  and  Griem,  1968, have 

already  calculated  the  lines  H    and  H     .    These  calculations  have  been 

o  7 

extended  to  H^^^  by  Bengtson,  Kepple  and  Tannich,  1969,  using  identically 

the  same  computer  program.    The  results  are  plotted  in  Fig,  15  and 

comparing  the  line  shape  for  the  higher  series  members,  for  which 

lower  state  interaction  becomes  negligible,  with  our  results  in  Fig.  13 

one  realizes  a  significant  difference  .    In  particular,  their  calculations 

-5/2 

do  not  reveal  the  AX  decay  in  the  near  line  wing  for  intensities 

smaller  than  about  l/lO  of  the  maximum  intensity  at  /\uu  =  0  which  is 
discussed  above.    It  should  be  pointed  out  that  the  ion  field  dependent 
cutoff,  which  has  been  introduced  by  Kepple  and  Griem,  1968,  to  account 
for  the  usually  neglected  exponential  in  Eq.  (VI.  4)  cannot  be  responsible 
for  it.    This  has  been  tested  in  our  calculations.    One  can  understand 
this  by  realizing  that  for  the  higher  series  members  the  effect  of 
dynamic  broadening  due  to  the  electrons  as  described  roughly  by  the 
constants  p^  in  Eq.  (X.  5)  turns  out  to  be  much  smaller  than  the  half- 
width  of  the  total  line,  which  is  essentially  determined  by  quasistatic 
broadening.  - 

As  another  interesting  result,  Fig.  16  shows  a  plot  of  a  calculated 

Lyman-p  profile  for  two  different  values  of  the  constant  B(B  =  1,  27  and 

3     2  , 

B  for  p     .    =  ^  +  — —  n  a  ),  which  allows  also  some  qualitative  state- 
mm  2  o 

ments  concerning  H^.    We  realize  that  changing  B  affects  the  very  line 
center,  where  the  profile  shows  the  two  humps  and  the  near  line  wing, 
but  it  does  not  change  the  intensity  around  the  halfwidth  significantly, 
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which  may  be  understood  as  an  effect  of  the  normalization.    This  is  in 
agreement  with  experimental  observations  of  Wende,  1967,  which  show 
that  the  calculations  of  Griem,  Kolb  and  Shen,  1962  overestimate  the 
near  line  wing.    It  also  explains  the  good  agreement  of  experimental  and 
theoretical  halfwidths  in  high  density  plasmas  (see  Gerardo  and  Hill, 
1966)  because  the  line  intensity  around  the  half  width  is  rather  insensitive 
to  the  exact  value  of  B. 

Finally  in  Fig.  17  to  19,  we  compare  the  unified  theory  cal- 
culations (solid  curves)  with  calculations  based  on  the  one-electron 
theory  in  order  to  see  how  far  into  the  line  center  the  asymptotic 
wing  expansions  as  given  in  Eq.  (XL  10)  or  (XI.  12)  may  be  used.  In 
all  Figures  the  short  vertical  line  indicates  the  position  of  the  outer- 
most, unperturbed  Stark  component  for  an  average  ion  field  B^^. 

which  is  given  by  Ati)  =  6      Auu .  (n,  q  =  n- 1 )  where  Auu .  is  defined  in 

av      1  1 

Eq.  (V.  14).    The  dashed  lines  correspond  to  the  one-electron  theory 

calculations  for  B  =  0  according  to  Eq.  (XL  12),  while  the  dash-dotted 

lines  give  the  results  for  6=R      according   to  Eq.  (XL  10),    First  of 

av 

all  we  realize  that,  as  expected,  the  one-electron  result  for  R  =  B 

av 

diverges  when  Auu  approaches  6     ^ti)  .{n,  q  =  n-1).    However,  in  all  three 

3.  V  X 

cases  we  see  that  for  frequencies 


Auu  >  5  B      Auu.  (n,  q  =  n-1)  (XILIO) 
av  1 

the  one-electron  theory  calculations  according  to  Eq.  (XL  10)  coincide 
with  the  unified  theory  calculations  to  within  1%  and  better.    For  slightly 
released  accuracy  requirements  one  may  also  apply  the  simpler  one- 
electron  theory  calculations  based  on  Eq.  (XI.  12)  with  3=0.  In 
particular  we  see  that  for  small  principal  quantum  numbers  the  useful 
range  is  very  much  larger  than  for  the  one-electron  theory  calculations 
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with  B  =         because,  for  p  =  0,  the  one-electron  theory  diverges  only 
at  AtA)  =  0.    We  also  realize  that  for  the  line  intensity  range  of  practical 
interest  both  asymptotic  wing  expressions  with  B  =  0  and  S  =  B 

av 

become  less  useful  with  increasing  principal  quantum  number. 
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APPENDIX  A 

PROGRAM  FOR  CALCULATING  THE  THERMAL  AVERAGE  F(t) 

This  section  gives  a  complete  listing  of  the  program  which  was 
used  to  calculate  the  thermal  average  F(t)  as  discussed  in  Sec.  IX.  In 
I  order  to  understand  the  program  the  following  explanations  may  be  of 
help. 

1.    Calculation  of   $    (s,  y,  x,  u) 

Function  PHI(  KAY)  calculates  $  ,  $      and  ?      as  defined  by 
the  Eqs.  (IX,  17)-(IX.  21).    In  order  to  assure  sufficient  numerical 
accuracy  for  computing  $    (s,  y,  x,  u),  series  expansions  have  been 
applied  whenever  one  of  the  different  g    (s,  y,  x,  u)  becomes  very  small. 
The  following  expansions  have  been  used  abbreviating 


w 


=    Yp    +  V  t^   /D  =  X    +  y    (1-x  ) 

=  p/  "Yp    +  v  t^  =  x/r 

/  2  2  '  ^  / 

=  vt  /  +  V  t      =    ^jl-x  y/r 


(A.l) 


,J"2       2  2  '  / 

and  Y  -  '^^/  TP    +vt  =us/ 


(a)    g^  {s,  y,  X,  u) 


a  <  0.  01    and  y  >  0 

2  3 

=  ^s   ^  1  +    V^^+v)      ^  +    y(16+80y+40y  +7y  )  ^4 


r^(l+Y)       (  8(1+y)^  128tH-Y)^ 
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(96+448y+1024y^+1520y^+880y^+252y^+33y^)  6 
1024(1  +  y) 


(A.  2) 


a  <  Oo  01  and  y  <  -  2 


2     (         ^   2       y^(16-5y)     4     y^32-192y  +  112y^-21y^) 

"""^    (         8(1-y)  128(1-y)  1024  (1-y) 

(A.  3) 


Y  <  0,  01 


-2  7  3 


1  -  w-Y  -  (3-llw  )  \  +  (9-17w  ) 


+  (31-350w^+447w^)  +  }  (A.  4) 

Y  >  100 


8 

4 


128     ■  I 


V2(l-wy  ^   L       (1+w)  r       5w-l  (3+4w-15w^) 
xu       •  I     "      2y     [  "      4y     "  3^2 


(ll-83w-59w^+195w^) 

3 

64y 


(b)   g^is.y,  X,  u) 


2 

(x/y)    <  0.  0002    and  y  >  0 


]! 


(A.  5) 
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3  6 

-  (187-443y+297y^-33y^  •  +  .  .  . 


128 

(A.  6) 

(x/y)^<  0.  0002    and  y  <  0 

.2  ,    V    2  .  .2  ,    V  4 


2 


2  xu 

(A.  7) 


I 1-yl  <  3. lO"^ 


! 


2  3 
+  (l-x^)    (l--^x^)(l-y)    +«..>  (A.  8) 


2.    Calculation  of  G(s,  x,  u) 

Function  GN(AX,  AXS)  calculates  G(s,  x,  u)  as  defined  by  the 
Eqs.  {IX.  22)  and  (IX.  24).     This  integral  has  been  solved  by  rewriting 
it  in  the  following  way. 


G(s,  X,  u)  =  G^(s,  X,  u)  +  G^Cs,  X,  u)  +  G^(s,  x,  u)  (A,  9) 


The  new  integrals    G  (s,  x,  u)  are  only  functions  of  $  (s,  y,  x,  u)  and 
turn  out  to  be 

1-R  1-R 

G,  =  U(x>  X  )  U(2  >  R)    I    $,  dy  +  U(x  >  x)  U(l  >  R+P)  |  $,dy 
1  o  J       1  o  J  1 

R  P 
"2  (A.  10a) 
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=  U{l-R  >  K) 


dy  with   K  = 


P  if   X  <  X 


R 

—  if  X  >  X 
2  o 


+1 


G2=U(x>x^)  j^U(2  >R)  j  ^^dy  +  \J(R>  2) J  ^^j 

1-R       1  -1  1 

+  U{x^>x)    j^U(l  >  R+P)J*  $^dy  +  U{R+P  >l)y  $2^^^ 


1-R 


1 

I 


G    =  I      $  dy  with  B  =  MAX  (1-R,  K)   and  K= 
B 


P 

- 1  if  X  >x 

o 

P  if  X  <  x 


and 
G 


=  U,x>x^).u(f>l)  -(f-l) 


4 


The  first  integral  is  split  up  at  most  into  three  parts 
-X  +x  1-R 

G  r 


lj=J     ?jdy+J     $jdy+J  l^dy 


K 


-X 


(A. 10b) 


(A.  11a) 


(A. lib) 


(A.  12) 


(A. 13) 


if  every  upper  limit  is  larger  than  the  lower  limit  using  a  different 
convenient  change  of  variables  in  every  part.    The  integration  is  per- 
formed in  all  cases  by  means  of  Weddle's  rule  (subroutine  WEDDLE). 
The  number  of  points  besides  the  initial  point  is  firstly  taken  to  be  6 
and  is  doubled  in  every  iteration  by  calculating  only  the  values  in  between 
the  old  values  until  the  integral  changes  less  than  a  preset  relative 
value  called  DG. 
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3.    Function  GUS  (Au)  calculates  the  x-integral  and  is  defined  by 
GUS(AU)=     =4=    u^»e"^    H(s,u)  (A.  14) 


where  1 

H(s,  u)  =    zfx-Vl-x     G(s,x,u)dx  (A.  15) 


I J  x^l-x^  '  ( 


o 


The  integral  is  performed  in  two  parts 

x^  -^n(x^) 

H(s,u)  =  /* ^1-z  'G{s,yJT^u)  dz  +  zf  e'^^-^l-e'^^  G(s,  e"^,  u)  dz 

o  o  (A. 16) 

by  means  of  Weddle's  rule  again  using  a  technique  as  described  before 
for  G(s,  X,  u).    The  new  test  parameter,  which  determines  the  number 
of  iterations,    is  called  ACC  and  is  given  in  the  main  program.  The 
old  test  parameter  DG  has  been  made  a  function  of  ACC,  x  and  u  in 
order  to  calculate  those  values  of  G(s,  x,  u)  with  the  highest  accuracy 
which  give  the  largest  contribution  in  the  final  velocity  integral.  For 
very  small  values  of  u*  s    the  integral  has  been  approximated  by 


Z  2 
H  (s,  u)  =  -  C  •  s 
a 


(A.  17) 


This  limit  does  not  lead  to  the  static  wing,  as  F(t)^    ^  in  Eq.  (IX.  28) 

does  because  of  the  lower  cutoff  x    .    As  explained  in  Sec.  IX.  this  cut- 

o 

off  is  necessary  in  order  to  avoid  rapid  fluctuations  in  the  integrals 

which  require  more  integration  steps  and  longer  computer  time.  For 

practical  calculations  x     is  chosen  small  enough  so  that  over  at  least 

o 

one  order  of  magnitude  in  t   the  calculated  F(t)  comes  as  close  as 

required  to  the  asymptotic  limit  F(t)  . 

t  -*0 
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For   us  >  2   the  following  relation  is  applied 

1 

2  f  xyi-x^'G^(s,x,u)dx=  -C  sf   -^iiLi|l£.)  dz 
J  J  1+z 


(A.  18) 


where 


+  ^(l-b  K  (b)  +  b^K  (b)  ) 
3  1  o 


u 


For  simplicity  x    has  been  set  to  zero  for  this  relation  without 
o 

affecting  the  final  result  noticably.    For  very  large    s    Eq.  (A.  18) 

leads  to  Eq.  (IX.  31).    However,  Eq.  (A.  18)  does  not  require    C   to  be 

small  as  Eq.  (  IX.  31)  does.     The  Raabe  integral  is  calculated  by  the 

function  SNZ  and  the  modified  Bess  el  functions  K    and  K    by  the 

o  1 

subroutine  BESMOD. 

4.    The  final  thermal  average  F(s,  n^^,  n^,  T)  as  defined  in  Eq. 
(IX.  23)  is  calculated  in  the  main  program  FSTEST.    The  best  results 
for  the  velocity  integral  have  been  obtained  by  Gauss's  quadrature 
formula  (function  GLQUAD).   The  values  FS  in  the  program  are  given 
by 


FS=   F(s,n,,  n  ,   T)   .  (A.  19) 

2nn  D-^  ^  ® 

e 

The  main  program  reads  in  the  temperature  T,  the  electron  density 

n  ,  the  quantum  number  n  ,  two  cutoff  parameters,  which  specify 
e  k 

the  lower  and  upper  cutoff  of  the  x- integral  by 

D      .    =  STRONG  •  p 

^"^^^  °  (A.  20) 

p  =  CUT  .  D 

max 
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where        is  given  by  Eq.  (IX.  5).    Finally  the  main  program  reads  in 
the  initial  value  of    s    for  computing  the  thermal  average  and  the  number 
of  values  which  proceed  according  to  s        =  s  •   10  .     The  program  also 

x^"r  J,  xC 

gives  the  asymptotic  thermal  average  leading  to  the  static  wing  which 
is  called  GS  and  calculates  the  relative  thermal  average  in  units  of 
this  asymptotic  value.    Furthermore,  the  function  KLOCK  provides 
a  means  to  test  the  computer  time  for  every  individual  value  of  FS. 
The  results  are  shown  and  discussed  in  Sec.  IX. 
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PROGRAM  FSTEST 
C  CALCULATION  OF   THE   THERMAL   AVERAGE  FS 

C  INTEGRAL   FS  OVER  U    IS   TESTED  FOR   DIFFERENT  L    IN  GAUSS  QUADRATURE 

COMMON/ PDS/ S/ PDCON /C 1 ♦CZ  » CON ♦ BCON /PDR AD/ BR AD / PDACC/ ACC ♦ DG » NGN 
COMMON /PDSTR/ STRONG 
EXTERNAL  GUS 

500   READ(60»100)    T EMP » DENS T Y , QNUM ♦ Q , CUT ♦ STRONG »    SS ♦  K 
IF ( DENSTY.EQ.O. )    CALL  EXIT 

WRI TE { 61  ♦  104 )    TEMP*   DENSTY*   QNUM,   Q,   CUT»  STRONG 

T   =  TEMP*l.E-4 

RELDEN   =   DENSTY^l .E-18 

CI    =   1.5   *  Q  *  QNUM 

C2   =   0.03043    ^   SORT ( RELDEN )/( T   *  CUT) 

WRITE    (61,105)  C1,C2 

NK   =  QNUM  *   0  +  0.01 

CON   =   Cl   *  C2 

BCON   =    1.414213562   *  CON 

CT3   =   CUT  3 

BRAD   =  7.67E-3*SQRT  (RELDEN/T  )*QNUM^-*2*STR0NG/CUT 

NUM  =  0 
ACC   =  3.E-4 
7   DO   1    I    =  1,K 
NUM  =   NUM  +  1 
SS   =   SS  *  10.0 
S   =   SS   /  CUT 
FSOLD   =  0.0 

GS   =   -1.671085516  ^    (CON*CUT*SS)  1.5 

GSAS   =   -(  CON->^SS  )  **2*    (  1  .  1 2  8  3  79 1 6  7 1  /  (  C2 *ST  RONG  )    -   1.)    *  CUT 
PRINT   200, SS,   GS,  GSAS 
DO  20   L   =  3*5 
LLL   =  KLOCK(O) 

FS   =  GLQUAD    ( GUS , 0 . 0 , 5 . 0 , L )    *  CT3 
FFGG   =  FS/GS 

DACC   =   ABS( (FSOLD  -  FS)/FS) 
LLL   =   KLOCK(O)    -  LLL 
PRINT   300,      FS,   FFGG,   DACC  LLL 
IF ( DACC.LT. ACC )    GO  TO  50 
20   FSOLD   =  FS 

50   PUNCH  400, SS,   FS,   DACC,   NK ,   CUT,   STRONG,    DENSTY,   TEMP,  NUM 
1  CONTINUE 
GO   TO  500 

100   FORMAT(    2E10.2,   4F10.1,    ElO.2,  15) 

104  FORMAT    (1H1»*   TEMPERATURE   =   *E14.5,    lOX,   *  DENSITY   =  *E14.5// 

1  *  QUANTUMNUMBERS  N  =   *F6.1,10X,*      (Nl   -  N2  )    =  ■»-F6.1// 

2  *   DEBYE  CUTOFF   FACTOR   =   *F8.3,10X,*   STRONG  COLLISION   FACTOR  = 

3  F8.3) 

105  FORMAT      (//*  CI   =   *   E14.6    ,    lOX       C2   =   *  E14.6) 

200   FORMAT    {///^  S   =   *E15.7,9X,*  GS  =   *E17.9,9X,*   GSAS   =  *E17.9/) 
300   FORMAT    (*   FS   =  ^E17.9,*   FS/GS   =   *E17.9,*   DACC   =  *E17.9, 

1        *  LLL   =  *I8/ ) 
400   FORMAT    (ElO.3,    E15.7,   ElO.2,    15,    F7.1,    F7.3,    2E11.3,  14) 

END 

C 

c 
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FUNCTION  GUS(AU) 

COMMON/PDS/S/PDCON/Cl »C2  jCON  » BCON/PDU/U/PDXM/XMl yXfAZ 
COMMON/ PDRAD/BRAD/PDACC/ AC C  »DG» NGN/ PDSTR/ STRONG 
DIMENSION   F( 192 )  ♦   H ( 192  ) 
GUS   =  0.0 
U   =  AU 

IF(S.EQ.O.O)  RETURN 

XMl   =  C2*STR0NG/U 

IF ( XMl .GE . 1 . )    XMl   =   1.    -  l.E-9 

XM2   =   XMl   +  BRAD 

IF(XM2.GE.l. )      XM2   =    1.    -  l.E-9 

CONS   =   CON  *  S 
U2   =   U  *  U 

FU   =   2.2567583342   *  U2   *  EXPF(-U2) 
US   =   U  *  S 

IF( (US/STRONG) .GT. 2. E-6)    GO  TO  5 
GUS   =   -CONS  *  CONS  *    (l./XM]    -   1 . )    *  FU 
RETURN 
5   NGN   =  0 
GD   =  0.0 

IF(US.LT.2. )    GO   TO  8 

NGN   =  1 

PA   =   2.   *  CON/U 

XMM   =   XM2   *  XM2 

SXM  =   SORT ( 1 .   -  XMM ) 

CALL   BESMOD    (PA,   F I  0 , F 1 1 , FKO , FK 1 ) 

GD  =   -CONS  *   SNZ(PA)    +  2 . * ( 1 . +PA* ( PA*FK0  -  FKl))/3. 
8   GUOLD   =  0.0 
GIOLD  =  0.0 
G20LD   =  0.0 
N  =  3 

DO   100  K   =  1»6 
GUS   =  GD 
DG   =  ACC/FU 
N   =   2   *  N 
NN   =  2 

IF(K.EQ.l )      NN   =  1 

AN   =  N 
ANN   =  NN 

I  F ( K. EG. 1 )    GO  TO  10 

DACC   =   ABS{ DGl /GUOLD ) 

I F ( DACC.LT. ACC )    GO  TO  40 

GO   TO  20 
10   FO   =  GN(0.0»1.0) 

DQl    =   XMl   *  XMl 
20   DQO  =  DOl/AN 

DQ  =   DQO  *  ANN 

Q  =   DOO  -  DQ 

DO   30   J  =  1»N,NN 

Q  =  Q  +  DQ 

X   =  SORT(Q) 

XS   =   SORT (1.   -  Q) 
30   F(J)    =   GN(X,XS)    *  XS 
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CALL   WEDDLE    ( DQO , N ♦ F » GUSl , FO ) 
40   GUS   =   GUS  +  GUSl 

IF(K.GE.6  )    GO  TO  50 

DGl   =   GUSl   -  GIOLD 

GIOLD   =  GUSl 

IF{ DGl. EG. 0. )    GO   TO  50 

DO  45   J   =  1»N 

L   =   N  +   1   -  J 
45   F(2   *   L)    =  F(L) 

50    I F ( XMl .GE .0.999999999 )    GO   TO  90 

IF(K.EQ.l )    GO  TO  60 

DACC   =   ABS{ DG2/GU0LD) 

IF(DACC.LT.ACC)    GO  TO  80 

GO   TO  6  5 
60   XO   =   SQRTd.   -  DQl) 

HO   =   GN(XM1»X0)    *   XO  DOl 

DYl    =  -LOGF ( XMl ) 
65   DYO   =  DYl/AN 

DY   =   DYO  ^  ANN 

Y  =   -DYl   +   DYO   -  DY 
Nl    =   N  -  1 

DO  70   J   =  1»N1»NN 

Y  =   Y  +  DY 
X   =  EXPF(Y) 
X2    =   X   *  X 

XS   =   SQRTd.   -  X2) 
DG   =  ACC/XS 
70   H(J)    =   GN(X»XS)    *   XS  *  X2 
H ( N )    =  0.0 

CA'.L   WEDDLE    (  DYO  » N  ♦  H  ♦  OUT  ♦  HO  ) 

GUS2   =   2.0   *  OUT 
80   GUS   =   GUS  +  GUS2 

IF(K.GE.6)    GO  TO  90 

DG2   =   GUS2   -  G20LD 

G20LD   =  GUS2 

IF(DG2.EQ.O. )    GO  TO  90 

DO  85   J   =  1»N 

L   =   N  +   1   -  J 
85   H(2   *   L)    =  H(L) 

90   DACC   =  ABS((GUS  -  GUOLD)/GUS) 

IF( DACC.LT.ACC  )    GO  TO  120 
100   GUOLD   =  GUS 
120   GUS   =   GUS  *  FU 

RETURN 

END 

C 
C 
C 

FUNCTION  GN{AX»AXS) 

COMMON/ PDS/S/PDU/U/PDCON/ CI ♦ C 2 ♦ CON , RCON/ PDXD/X , XS ♦ X 2  » DX 2 
C0MM0N/PDACC/ACC»DG»NGN/PDXM/XM1 ♦XM2 
DIMENSION     F(768),   G{768)»   H(768)»  D(768) 
GN   =  0.0 
X   =  AX 
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XS   =  AXS 
X2    =   X   *  X 
0X2   =    1.    -  X2 
GiMOLD   =  0.0 
GIOLD  =  0.0 
G20LD  =  0.0 
G30LD  =  0.0 
G40LD   =  0.0 
R   =   U   *   S   /  XS 
Rl   =   1.   -  R 
R2   =   0.5   *  R 
PI   =  -1. 
P2   =  -1. 
N  =  3 

DO   100  K   =  1*8 
GN   =  0.0 

N   =    2   *  N 
NN   =  2 

IF(K.EQ.l )    NN   =  1 

AN   =  N 
ANN   =  NN 
QA   =  -R2 

IFCX.GE.XMl)    GO   TO  5 
IF(K.EQ.l)    PI   =  SQRT(XM1*XM1 
QA   =  PI 
5    IF(Rl.LE.Pl)    GO   TO  50 
QB   =  -X 

IF(Rl.LT.-X)    OB   =  Rl 
IF(QA.GE.QB)    GO  TO  20 
IF(K.EQ.l )    GO  TO  7 

DGG   =  ABS(DG1/GN0LD) 
IF( DGG.LT.DG)    GO  TO  13 

GO  TO  8 

7  DQl   =    (QB  +  R2)  ( 1./3.  ) 

8  DQO  =  DQl/AN 
DQ   =   DQO   *  ANN 
0  =   DQO  -  DQ 

DO   10   J   =  1»N»NN 

Q   =   Q  +  DO 

02   =  Q  *  0 

Y   =   QR   -   0  *  02 
10   F  (  J )    =   PHI (  1 »Y)    *  Q2 

CALL  WEDDLE    ( DQO » N » F » OUT  ,  0  .  0 

Gl   =   3.0  OUT 
13   GN   =  Gl 

IF(K.GE.8)    GO  TO  18 

DGl   =   Gl   -  GIOLD 

GIOLD  =  Gl 

IF(DG1.EO.O» )    GO  TO  18 

DO   15   J   =  1»N 

L   =   N  +   1   -  J 
15   F(2   *   L)    =  F(L) 
18   QA   =  -X 
20  QB   =  X 


IF(R1.LE«X)    QB   =  Rl 

IF(QA.GE.OB)    GO   TO  35 

IF(K.EQ.l )    GO  TO  23 

DGG   =  ABS(DG2/GN0LD) 

IF(DGG.LT.DG)    GO   TO  31 

GO   TO  25 
23   FO   =   PHI ( 1 ,QA ) 

DY2   =   QB  -  OA 
25   DYO   =  DY2/AN 

DY   =   DYO   *  ANN 

Y  =   OA   +  DYO  -  DY 
DO   30   J   =  1»N»NN 

Y  =   Y   +  DY 

30  G ( J )    =   PHI { 1 ♦ Y ) 

CA'-L  WEDDLE    (  DYO  »  N  ,  G  »  G2  »  FO  ) 

31  GN   =  GN  +  G2 
IFIK.GE.S)    GO  TO  33 
DG2   =   G2   -  G20LD 
G20LD   =  G2 

IF(DG2.EO.O. )    GO  TO  33 

DO  32  J  =  1»N 
L   =   N  +   1   -  J 

32  G(2   *   L)    =   G( L) 

33  OA   =  X 

35  IF(QA.GE.Rl)    GO   TO  50 
IF(K.EQ.l)    GO  TO  36 

DGG   =  ABS(DG3/GN0LD) 
IF(DGG.LT.DG)    GO   TO  42 
GO  TO  37 

36  HO   =   PH I  (  1  ♦Rl  )    *   Rl   *  Rl 
QO   =  l./Rl 

DQ3   =    l./QA  -  QO 

37  DQO   =  DQ3/AN 
DO   =   DQO  *  ANN 

Q  =   QO  +  DQO  -  DO 
DO  40   J   =  1»N»NN 
0  =   Q  +  DQ 

Y  =  l./Q 

40   H( J)    =   PHI ( 1 »Y)    ^   Y   *  Y 

CALL   WEDDLE    ( DQO ♦ N ♦ H ♦ G3 ♦ HO ) 
42   GN   =  GN  +  G3 

IF(K.GE.a)    GO  TO  50 

DG3   =   G3   -  G30LD 

G30LD  =  G3 

IF(DG3.EQ.O. )    GO  TO  50 

DO  45   J   =  1»N 

L   =   N  +   1   -  J 
45   H(2   *   L)    =  H(L) 
50    IF(R.LT,l.E-6)    GO   TO  80 

OA   =  Rl 

IF(X.GE.XM2)    GO  TO  55 

IFCK.EQ.l)    P2   =   SORT(XM2*XM2   -   X2 ) /XS 
55    IF(R1.LT,P2)    QA   =  P2 
IF(K.EQ.l)    GO  TO  57 
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DGG   =   ABS ( DG4/GN0LD ) 
IF( DGG.LT.DG)    GO  TO  63 
GO  TO  58 

57  DA   =  PHI(2»QA) 
DY4   =    1.    -  QA 

58  DYO   =  DY4/AN 
DY   =   DYO  *  ANN 

Y  =   QA   +   DYO  -  DY 
DO   60   J   =  1»N»NN 

Y  =   Y   +  DY 

60   D{ J)    =   PHI (2»Y) 

CALL   WEDDLE    ( DYO » N ♦ D » G4 ♦ DA ) 
63   GN   =   GN   +  G4 

IF(K.GE.8)    GO  TO  70 

DG4   =   G4   -  G40LD 

G40LD   =  G4 

IF ( DG4.EQ.0. )    GO  TO  70 

DO  65   J   =  1»N 

L   =   N  +   1   -  J 
65   D(2   *   L)    =  D(L) 
70    IF(R2.LE.l. )    GO   TO  80 

IF(K.EQ.l)    G5   =    ( R2-1. )*PHI (4»Y) 

GN   =  GN  +  G5 
80   DGG  =   ABS((GN  -  GNOLD)/GN) 

IF( DGG.LT.DG)    GO  TO  120 
100   GNOLD   =  GN 

120    IF(NGN.EQ.l)    GN   =   GN  -  G5 

RETURN 
END 


FUNCT ION   PHI ( K  »AY ) 

COMMON/ PdS/S/PdU/U/PDCON/ CI » C 2  » CON , BCON/ PDXD/ X  » XS , X 2  » 
Y   =  AY 
PHI    =  0.0 
FAC   =  0.0 

IF(K.GT.2)    GO  TO  30 

RI    =   l./SQRT(X2   +   Y  *  Y  -i^   DX2  ) 

IF{K.GT.l )    GO  TO  10 

C   =   X   *  RI 

W  =   XS  *   Y  *  RI 

G  =   u  *   S  *  RI 

IFCC.GT.O.Ol )    GO  TO  3 

IF ( Y.LT.O . )      G   =  -G 

AG   =   1  .  /ABS ( 1 .   +  G) 

IF( AG.GT.2. )    GO   TO  70 

GI    =    (C  *   AG)  2 

IF(G.LT.-1.)    GO   TO  1 

FAC   =   CON   *   S  *   RI   *  RI   *  AG 

ARG   =   FAC   *    (1.    +   GI    *   G   *   0.125   *    (4.    +  G)) 
I  F ( C.LT. 1 .E-3 )    GO   TO  40 

AG   =( ( (0.0546875*G+0.3125)*G+0,625)*G+0.125)*G*GI*GI 

GO   TO  35 
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1   FAC   =   2.0   *  CON/{X*U) 

ARG   =   FAC  *    (1.   -  GI    *  0.125   *  G  *  G) 
IF( C.LT. 1 .E-3 )    GO  TO  40 

AG   =   -(0.  125  +  0. 0390625*G)    *  GI    *  GI   *  G-*^*3 
GO  TO  35 
3    IF(G.GT.0.01 )    GO  TO  5 
FAC   =   CON   *   RI    *  RI    *  S 
W2   =   W  *•  W 

ARG   =   FAC  *    (({1.375*W2   -   0.375)*G  -  W)*G  +  1.) 
IF ( G.LT.^ .E-4 )    GO   TO  40 

AG   =    { 3. 4921875*W2-2. 734375 ) *W2+0. 2421875 
AG   =    (AG*G  -    (2.12  5*W2   -   1.125  )-^W)        G  **  3 
GO   TO  35 
5    IFCG.LT.  100.  )    GO   TO  8 

FAC   =   BCON  *   SQRT(1.   -  W)/(X*U) 
GI    =  l./G 

ARG  =   FAC  *    (  (  (  0.625-!^W-0.  125  )  ^GI-0.5  )  «-GI*  (  1  .+W  ) +  1  .  ) 
IF(GI  .LT.3.E-4)    GO  TO  40 

AG= ( (1. 52 34375* W- 0.4609375) *W-0.648  4  3  75 ) *W+0. 08593 75 
AG  =  { AG*GI+ (( 0. 25-0. 9375*W ) *W+0. 1 875 ))"( 1 .+W ) *GI**3 
GO  TO  35 

8   A  =    (1.   +  G*W)/SQRT(1.   +  G-^(2.*W  +  G)) 
GO   TO  20 
10   AG   =  ABStY) 

BY   =   ABS( 1.    -  Y) 

IF(AG.LT.l.E-lO)    GO   TO  15 

GI    =    (X/Y)  2 

IF ( GI .GT.2.E-4)    GO   TO  15 

IF{Y.LT.O.)    GO   TO  13 

FAC   =   CON   *   BY   /    (AG  *  U) 

ARG  =   FAC  *    ( l.-0.125*( l.+Y)*(3.-Y)*GI ) 
IF(GI.LT.l.E-6)    GO  TO  40 

AG   =    (  (0.0546875*Y-0. 328125  )*Y  +  0. 2421875  )  *(GI*(  l.+Y)  )  «-*2 
GO  TO  35 
13   FAC   =   2.0   *  CON/(X*U) 

ARG   =   FAC  *    ( 1  .-0. 125*GI*BY**2 ) 
IF(GI .LT.l.E-6 )    GO  TO  40 

AG   =    (( 2. -Y ) *0. 0390625*Y+0. 0859375 )*( GI*BY ) **2 
GO   TO   3  5 
15    IF(BY.GT.3.E-3)    GO  TO  18 
I F ( BY.LT. 1  .E-11  )  RETURN 
FAC   =   CON  *   XS  *   BY   /  U 

ARG   =   FAC   *    (  l.  +  (  (  1  .-1  .  375^^X2  )*BY   +  l.)*DX2*RY) 

IF(BY.LT.2.E-4)    GO   TO  40 

AG   =    (1.   -  2.125*X2)*BY*(BY*DX2)  2 

GO  TO  35 
18A=    (X2+Y*   DX2 )    *   R I 
20    IFCA.GE.l.)  RETURN 

ARG   =   BCON  *   SQRTd.   -  A)/(X*U) 

GO   TO  40 
30    IF(X.FQ.O.)  RETURN 

ARG   =   CON   *   2.0   «•   XS   /    (X   ^-  U) 

GO   TO  40 
35   ARG   =   ARG   +   FAC   *  AG 
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40    IF( ARG.LT.O. )    GO  TO  60 

IF(ARG.GT.   0.05)    GO   TO   5  0 

PHI    =   -(0.5   -   ARG*ARG/24. ) *ARG*ARG 

RETURN 
50   PHI    =   COS(ARG)    -  1.0 

RETURN 

60   WRITE    (61»1^U)    ARG»    FAC»    K,    X»   U,    Y,  S 
100   FORMAT    (^^   ERROR      NEGATIVE   VALUE   OF   ARG    IN   PHI»   ARG   =  *E17.9, 

1  *   FAC   =   *E17.9/*   K   =   *I1»*   X   =   *E17.9,*  U   =  *E17.9, 

2  *  Y   =   *E17.9,*   S   =  -^E17.9) 
RETURN 

70  WRITE    (61*200)  G 
200   l-ORMAT    (*   PHI    NOT  ACCURATE   ENOUGH    ♦   G  =  *E17.9) 
CALL  EXIT 
END 

C 
C 
C 

SUBROUTINE  WEDDLE    ( DX ♦   N,   F»   A,   FO ) 

C  F   IS   THE   FUNCTION   TO  BE    INTEGRATED  BY  WEDDLES  RULE 

C  FO    IS   THE   VALUE  OF   THE   FUNCTION   TO  BE    INTEGRATED  AT   SOME  STARTING 

C  POINT   WHICH    IS   NOT    INCLUDED    IN   THE    INPUT   ARRAY  F 

DIMENSION  F(N) 

1  A   =  0.0 

2  K   =   N   -  1 

3  DO   15    I    =  1»6 

4  SUM   =  0.0 

5  DO  6   J   =    I  »   K  »  6 

6  SUM   =   SUM   +  F(J) 

7  GO   TO    (8»    1^»    12»    10,    8,    14),  I 

8  A   =   A   +   5.0   *  SUM 

9  GO   TO  15 

10  A   =   A   +  SUM 

11  GO   TO  15 

12  A   =   A   +   6.0   *  SUM 

13  GO   TO  15 
14A=A+2.0*  SUM 

15  CONTINUE 

16  A   =   0.3   *   DX   *    ( A   +   FO   +  F(N)) 

17  RETURN 
END 

C 
C 

c 

FUNCTION  SNZ(X) 

C  SNZ(X)    CALCULATES   RAABES    INTEGRAL   OVFR   S I N ( X*Z ) / ( 1 . +Z*Z )  DZ 

DIMENSION   AA(6),    BB(4),    CC ( 4 ) 

DATA( ( AA( I ) , 1=    1,6)    =      -0.57721566,   0.99999193,  -0.24991055, 
1  0.05519968,   -0.00976004,  0.00107857) 

DATA ( ( BB( I ) ♦ I =    1,4)    =      8.5733287401,    18.0590169730,  8,6347608925, 
1  0.2677737343) 

lATA( (CC( I ) ♦ 1=    1,4)    =      9.5733223454,    2  5.6329561486, 
1        21.0996530827,  3.9584969228) 

SNZ   =  0. 
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IF(X.FO.O.)  RETURN 
A  =   ABS(X ) 

I  F ( A.GT.40. )    GO  TO  100 
A2   =  A   *  A 
EMA   =  EXPF(-A) 
IF(A,GT.0.2>    GO   TO  10 

EEA   =    ( ( { (A2/72.+1. )*A2/42.+l. )*A2*0.05+1. )*A2/6.+l. )*A 

GO   TO  20 
10   EEA   =   0.5   *    (EXPF(A)    -  EMA) 
20    IF( A.GT.l. )    GO  TO  30 

SUM   =   AA( 1 ) 
Z   =  A 

DO  25   J   =  2»6 
SUM   =   SUM  +  AA(J)    *  Z 
25   Z   =   Z   *  A 

FIT    =   SUM   -  LOGF(A) 
SNZ   =   FFA  *  FIT 
GO   TO  40 

30   SUM  =    (  (  ( A  +  BB ( 1  )  ) *A  +  BB ( 2 )  ) *A  +  BB { 3 )  ) *A  +  BB ( 4 ) 
SUMM   =    { ( ( A+CC ( 1 ) ) *A+CC ( 2 ) ) *A+CC ( 3 ) ) *A+CC (4 ) 
FIT   =    SUM/(SUMM  *  A) 

SNZ   =   0.5   *   FIT   <-    (1.   -   EXPF  ( -2.-!^A )  ) 
40   PROD  =  A 
SUM  =  A 

DO  45   J   =  1*200 
AJ   =   2   *  J  +  1 

PROD  =   PROD  *  A2/ ( AJ* ( AJ-1 . )  ) 
SUM  =   SUM   +  PROD/AJ 
PT   =   PROD  *  l.F+10 
IF(SUM.GT.PT )    GO  TO  50 
45  CONTINUE 

50   PT   =   SNZ  +  SUM  *  EMA 
SNZ   =  SIGNF(PT»X) 
RETURN 
100   PROD  =  l./X 

A2   =   PROD  *  PROD 
SUM  =  PROD 
DO   1 10   J   =   1  » 100 
AJ   =   2   *  J 

PROD  =   PROD  *  A2   *   AJ  *    (AJ  -  1.) 

SUM   =   SUM  +  PROD 

IF{AJ.GF.A)    GO  TO  150 

PT   =   ABS(PROD  *  l.E+10) 

IF( ABS( SUM) .GT.PT)    GO  TO  150 
110  CONTINUE 
150   SNZ   =  SUM 

RETURN 

END 

C 
C 

C 

SUBROUTINE  BESMOD ( X ♦ F I  0 ♦ F I  1 ♦ FKO » FK 1 ) 
C  BESMOD  CALCULATES   THE  MODIFIED  BESSELFUNC T I ONS    10,    II,   KO  AND  Kl 

C  BY  MEANS  OF  POLYNOMIAL   APPROXIMATIONS  AS  GIVEN   IN  THE 
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NBS  HANDBOOK  OF  MATHEMATICAL   FUNCTIONS*     PAGE  378 
DIMENSION   A(7)»B(9)»C(7),D(9)  ♦£{?), F(7)»G(7),H(7) 
DATA    ((A(I)»I    =    1,7)    =   0.0045813,    0.0360768,  0.2659732, 
1        1.2067492,   3.0899424,      3.5156229,  1.0) 
DATA    {(B(I)»I    =    1,9)    -   0.00392377,   -0.01647633,  +0.02635537, 

1  -0.02057706,   0.00916281,   -0.00157565,   0.00225319,  0.01328592, 

2  0.39894228) 

DATA    ({C(I)»I    =   1,7)    =   0.00032411,   0.00301532,  0.02658733, 
1        0.15084934,    0.51498869,    0.87890594,  0,5) 
DATA    ((D(I),I    =    1,9)    =   -0.00420059,    0.01787654*  -0.02895312, 

1  0.02282967,   -0.01031555,   0.00163801*  -0.00362018,-0.03988024, 

2  0.39894228) 

DATA    ((E(I),I    =   1,7)    =   7.4E-6,    1.075E-4,   0.00262698*  0.0348859* 
1        0.23069756*   0.4227842*  -0.57721566) 

DATA    ((F{I)*I    =    1*7)    =   0.00053208*   -0.00251540,  0.00587872, 
1       -0.01062446,   0.02189568,   -0.07832358,  1.25331414) 

DATA    ((G(I)*I    =   1*7)    =   -4.686E-5*   -0.00110404,  -0.01919402* 
1        -0.18156897*   -0.67278579*    0.15443144*  1.0) 

DATA    ((H(I)*I    =    1*7)    =   -0.00068245*   0.00325614*  -0.00780353, 
1        0.01504268,   -0.0365562*   0.23498619*  1.25331414) 

XS   =  0.0 

XE   =  0.0 

IF(X.LT.2. )    GO  TO  10 

XS   =  SORT(X) 
XE   =  EXPF(X) 
10    IF(X.GT.3.75 )    GO  TO  20 

Y  =  X/3.75 

Y2    =   Y   *  Y 

FIO   =    (  (  (  (  (  Y2*A ( 1 )+A ( 2  )  ) *Y2+A { 3  )  ) *Y2  +  A (4  )  ) *Y2  +  A( 5 )  ) * Y2  +  A ( 6  )  ) * Y2  +  1 . 

FIl   =(((((  Y2*C ( 1 ) +C ( 2  )  )*Y2  +  C(3)  ) *Y2+C ( 4 )  ) *Y2  +  C ( 5  )  ) *Y2+C ( 6  )  ) *Y2  +  0.  5 

FIl    =   FIl   *  X 

GO   TO  40 
20   Y   =  3.75/X 

FIO   =  8(1) 

FIl   =   D{ 1  ) 

DO   30   K  =  2*9 

FIO  =  FIO  *  Y  +  B(K) 
30   FIl   =   FIl   *   Y  +  D(K) 

XEX   =  XE/XS 

FIO   =   FIO  *  XEX 

FIl    =   FIl   *  XEX 
40    IF(X.GT.2.)    GO  TO  60 

Y  =  X/2. 
Y2    =   Y   *  Y 
FKO   =   E (  1  ) 
FKl   =   G(  1  ) 

DO   50   K   =  2*7 

FKO  =  FKO  *  Y2  +  E(K) 
50   FKl    =   FKl   *   Y2   +  G(K) 

XEX    =  LOGF(Y) 

FKO   =   FKO  -  XEX  *  FIO 

FKl    =   FKl/X   +   XEX   *  FIl 

RETURN 
60   Y  =  2.0/X 
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FKO   =   F( 1  ) 

FKl    =   H{ 1  ) 

DO  70   K   =  2»7 

FKO   =   FKO   ^   Y   +  F(K) 

FKl    =   FKl   *   Y  +  H(K) 

XEX   =    1 ./ { XS   *  XE) 

FKO   =   FKO  *  XEX. 

FKl   =   FKl   *  XEX 

RETURN 

END 


FUNCTION  GLQUAD(F,A,B,L) 
GAUSSIAN-LEGENDRE  QUADRATURE 
40  OR   80   NODES  FOR   L   =   1»  2» 
COMMON/GLQDAT/Xl  (2)  »W1(2)  »X2 
X5 ( 20  )  »W5 ( 20  )  ♦XS { 40  )  ♦we ( 4 


DATA 
DATA 
DATA 
DATA 
DATA 


(  (  XI  (  I 
(  (  Wl  (  I 
(  ( X2 (  I 
(  (  W2  (  I 
(  ( X3 (  I 


1  .973906528b 
DATA      (  ( W3 (  I 

1  ,0666713443 
DATA      ( ( X4{ I 

1  .9931285992 

2  .6360536807 
DATA      (  ( W4{  I 

1  .0176140071 

2  ,1181945320 

DATA      ( ( X5 ( I 

1  ,9982377097 

2  .9020988070 

3  .6719566846 

4  ,3419940908 
DATA      ( { W5 ( I 

1  ,0045212771 

2  ,0334601953 

3  ,0574397691 

4  ,0728865824 

DATA  ((X6(I 

1  ,9995538227 

2  ,9749091406 

3  ,9132631026 

4  .8169541387 

5  .6896376443 

6  .5361459209 

7  .3623047535 

8  .1747122918 
DATA  (IW6(I 

1  .0011449500 

2  .0086839453 


I  =  1^2)  = 
I  =  1 ^2  )  = 
I  =  1^3)  = 
I  =  1^3)  = 
I  =  1^5)  = 
.8650633667 
I  =  1^5)  = 
.  1494513491 

I    =  1»10) 
.9639719273 
.5108670020 

I   =    1 ♦ 10 ) 
.0406014298 
.1316886384 

I    =  1^20) 
.9907262387 
.8659595032 
.6125538897 
.2681521850 

I    =  1^20) 
.0104982845 
.0387821680 
.0613062425 
.0747231691 

I  =  1^40) 
.9976498644 
. 9654850890 
.8966755794 
.7938327175 
.6608598990 
. 5028041 1 19 
. 3256643707 
.1361640228 

I    =  1^40) 
.0026635336 
.0101617660 


OF  F   FROM  A 
3^   4^    5^  6 
(  3  )  ♦WP  (  3  )  »X3 

•0  ) 

8611363116 
3478548451 
9324695142 
1713244924 

6794095683 

2  190863625 

9122344283 
3737060887 

0626720483 
1420961093 

.9772599500 
8246122308 
5494671251 

.1926975807 

.0164210584 
.  0438709082 
.0648040135 
.0761 103619 

.9942275410 
9545907663 
.8787225677 
.7695024201 
6310757730 
.4686966152 
.2885280549 
0974083984 

,0041803131 
,0116241141 


TO   R   WITH  4 

( 5 )  ♦ W3 ( 5 )  ♦X4 

.3399810436 ) 
.6521451549) 
.6612093865 
.  3607615730 

.4333953941 

.2692667193 

.8391169718 
.227785851 ] 

.0832767416 
.  1491 729865 

.9579168192 
.7783056514 
.4830758017 
.1160840707 

.0222458492 
.0486958076 
.0679120458 
.0770398182 

.9892913025 
.9422427613 
.85943  14067 
.  7440002976 
,6003306228 
,4338753708 
.2509523584 
.0585044372 

.0056909225 
.0130687616 


♦6»10^20» 
( lO )  ♦W4(  10 ) 


.2386191861 \ 
.4679139346 

.1488743390 ! 

.2955242247 : 

.7463319065. 
.0  76  5265211: 

.1019301198. 
.1527533871 : 

.9328128083- 
,7273182552 . 
,4137792044. 
.0387724175 : 

,0279370070< 
.0532278470- 
.07061 16474, 
,0775059480 ; 

.9828485727, 
.9284598772 , 
.8388314736, 
.7173651854. 
.5686712681 , 
.3983934059, 
.2129945029, 
.0195113833 ] 

,0071929048, 
,0144935080, 


90 


3  .0158961836*  . 

4  .0225050902,  . 

5  .G28259816l»  . 

6  .0329419394,  . 

7  .0363737499,  , 

8  .038424993U,  , 
TO  =  ( A  +  B) /2. 
Tl  =  ( R  -  A) /2. 
Y  =  0. 

GO  TO  (1,2,3,4, 

1  DO   10   K   =  1»2 
10   Y=Y+W1 { K ) * ( F ( TO 

GO  TO  100 

2  DO  20  K  =  1»3 
20   Y=Y+W2 (K)*( F( TO 

GO  TO  100 

3  DO  30   K   =  1»5 
30   Y=Y+W3 ! K) * ( F ( TO 

GO  TO  100 

4  DO   40   K   =  1»10 
40   Y=v+W4 ( K ) * ( F ( TO 

GO   TO  100 

5  DO   50   K   =  1»20 
50   Y=Y+W5 ( K) * ( F ( TO 

GO   TO  100 

6  DO  60   K   =  1»40 
60   Y  =  Y  +  W6(K)*-(F(T0 

100  GLQUAD=Y*T1 
RETURN 
END 


0172746521 ♦ 
0237318829, 
0292883696  » 
0337332150, 
0368977146  » 
0386617598  » 


.0186268142  » 
.0249225358  , 
.0302723218  » 
.0344731205  » 
.0373654902  ♦ 
.0388396511  , 


5,6)  L 

-Tl*Xl(K) )+F(TO+Tl*Xl (K) 
-Tl*X2(K) )+F ( T0+T1*X2 (K) 
-T1*X3 ( K ) ) +F { T0+T1*X3 ( K ) 
-Tl*X4(K)  )+F(T0  +  Tl^^X4(K) 
-Tl*X5  (  K  )  )  +F  (  T0  +  T1^^X5  (  K  ) 
-Tl*X6(K) ) +F ( T0+T1*X6 ( K ) 


.0199506109 , 
.0260752358 , 
.0312101742 , 
.0351605290 , 
.0377763644, 
.0389583960 , 


.0212440261 , 
.0271882275 , 
.0321004987, 
.0357943940, 
.0381297113, 
.0390178137) 


i 
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APPENDIX  B 

THE  LARGE  TIME  LIMIT  OF  THE  THERMAL  AVERAGE  F(t) 

In  Eq.  (IX.  31)  the  large  time  limit  of  the  thermal  average  has 
been  given,  which  is  of  the  form 


This  form  has  been  obtained  by  most  modern  impact  theories.  The 
additive  constant  B  varies  depending  on  what  type  of  cutoff  has  been 
usedo    In  the  following  we  derive  the  different  constants  B  for  the 
different  cutoff  procedures  which  have  been  used  and  compare  them 
with  the  numbers  given  in  the  literature. 


F(t,  n  ,  n  ,   T),  differ  essentially  in  three  respects,  namely  by  the 
k  e 

upper  and  lower  limits  of  the  p -integral  and  by  the  limits  of  the  t'- 
integral  in  Eq.  (VIII.  4).    Based  on  the  completed  collision  assumption 
(Baranger,  1962),  the  limits  of  the  latter  integral  are  usually  extended 
from  -00  to  +  <».     This  approach,  however,  is  not  quite  consistent  with 
the  cutoff  at  the  Debye  length,  which  would  rather  require  the  integral 
to  go  from  -T  to  +  T  as  done  in  this  report  (T  is  defined  by  Eq.  (IX.  7)). 
We  therefore  have  to  investigate  the  following  integrals: 


F(t)^        =  At  (B-^n  (4C^)). 


(B.l) 


The  various  methods  to  evaluate  F  ,the  large  time  limit  of 


00 


a.) 


(case  a) 


00 


1 


X 


(B.2) 


mm 


92 


and 


b.  )   for/      PV^(t')dt'   ^/   PV^(t')dt'  (case  b) 


rj    PV^(t')  dt'   

1  du  ==,  u  e       I    dx  X     cos   I  J  - 


"~b  Z 
F    =  2^n  (a.  D)  V    t  /du  =^ 

00  e  av 


(B.3) 


X 

mm 

where 


X    .    =  p     .    /  (a-  D)  .  (B.  4) 

mm  mm 


The  factor  a  =  p        /D   is  usually  taken  to  be  one  and  has  in  some 
max 

papers  (e.  g.  ,  Griem  et  al.  1962)  been  varied  to  1.  1.  or  to  0.  606 
as  proposed  by  W.  -R.  Chappell,  J.  Cooper  and  E.  Smith,  1969. 
As  a  lower  cutoff  we  consider  in  particular  the  three  cases  of 

p     .    =  0,  p     .    -  K  and  p     .      =  3(nq-n'q')>^      2CD/u  by  setting 
mm  mm  mm        ^  >  i       7  & 

(B.5) 


mm      a  u 

In  the  following  we  will  set  a=l.    In  order  to  recover  the  dependence 
on  the  upper  cutoff  parameter  a,  we  only  have  to  replace  in  all  the 
following  relations  C  by  C/a. 

First  of  all,   one  realizes  that  with  x    .    <  1    the  lower  limit  on  the 

mm 

u- integral  is  given  by 

u   =  b-  2C  .  (B.  6) 

o 

Hence,  we  have  to  evaluate  the  following  two  integrals 

2 


=y     du  Je'^^^j     dxx  l^cos  ^  ^-x^'j-lj  (B.  7) 

o/  u 


,a 


u 

o 
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and 


dx  X 


u 


u  /u 


cos 


(B.  8) 


o  o- 

The  second  integral  can  be  simplified  after  a  change  of  variables  and 
a  partial  integration  to 


du 


(B.9) 


After  expanding  the  cosine  and  another  change  of  variables  we  have 

?. 

2       00  k 

^     r-n     /?r.\2k    /        ^  o 

dz  (B. 10) 


4     Z-^     (2k)!    U     1  / 
k=l  ^  o  /  A 


which  can  be  expressed  in  terms  of  exponential  integrals 


2k 


b  C 
I    =  -  E 


,2         ""o  (-1)  /2C\ 

k=2  ^    o  / 


E  (u  )   .  (B.ll) 

k  o 


With  the  lower  cutoff  parameters  stated  above  (b  <  1)  and  typical 

densities  and  temperatures  of  interest  one  usually  has  u    <  0.  1  (see 

2  2  ° 

Eq.  (IX.  21)).    Since  for  k>  2    E(u     )  =  l/(k-l )  +  0(u  )   one  obtains 

2  ^    ^  ° 

to  lowest  order  in  u 

o 


2k-2 


(B.  12) 


.2  r 


k-l 


1 


(2k)!     (2k-l)!  (2k-2)(2k-2)! 
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+  2 


fr)f""(f)(f)|-H(-^)-(t)|]- 


which  yields 


3{y-1)  +       (4C^)  +  2K  (  ^  +  0  (u")  |  (B.  13) 


where    K  is  defined  as 


,  \       1  -cos  z    ,      sin  z  ,  .  ,  .X 

K  (z)  =  —  +  Ci  (z)  (B.  14) 

L.  Z 

z 

and  Ci  is  the  cosine  integral.  Eq.  (B.  13)  was  obtained  already  by  Shen 
and  Cooper,   1969.     Their  constant  A  is  identical  with  our  constant  2C. 

The  other  integral  I    of  Eq.  (B.  7)  one  can  obtain  by  evaluating 

Al  =  I     du   u  e       I      dx  X      cos    I  '%/1-x  l-cos 

J  J         L     V""'  / 


so  that 


=      +  Al   .  (B.  16) 


If  we  again  expand  the  cosine  functions,  aI  can  be  given  by 

00  y,        k-1  .         no  _  1 

.     \->  (2C)        '^"^       k!  (-l)J  /    ,      3-2k  -u  /  l-2j 

\  V"  (B.17) 

— —  e 
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2k-2 


SI  /2c 
(2k)!  ^ 


k=2 


2C\ 


E       (u  )-E  ( 
k-1     o  k 


o      I  k-1 


dz 


2k-2  ^ 

4  2  k!      /  2C\  (-l)-^ 

""o  (2k)!     I  u    I  Z^j!(k-j)!| 

k.3  \  0/ 


(2j-2) 


k-1     o       k-j  o 
(B.17) 


which  gives  us  to  lowest  order  in 


u 


C 


1  +  0(u  ) 

o 


(B. 18) 


This  means  that  for  the  same  lower  cutoff  case  a.  )  and  b.  )  as  defined 
in  the  Eqs,  (B.  2)  and  (B.  3)  differ  only  by  a  constant  1  in  their  additive 
constants  B.     As  a  result  we  have 


F    =  -    (-l-(nq-n'q')-^l     n  t 
00  12  m  I  e 


B-;,n  (4C  )     (B.  19) 


'1 


where  the  constant  B  for  the  different  cutoff  parameters  is  compiled  in 
the  following  Table  2. 


Table  2.     The  constant  B  for  different  cutoff  parameters. 


+  T 

/     PV^(t')  dt' 
-  T 

+00 

/    PV  (f)  dt' 

J  " 

-00 

p     .    =  0 
mm 

0.  27 

1.27 

p     .    =  ^ 
mm 

0.  23  -  0.  27 

1.  23  1.27 

p     .     =  3n  ^ 
mm  k 

-1.  66 

-  0.  66 
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In  order  to  compare  our  results  with  the  numbers  given  in  the  lit- 
erature we  rewrite  Eq.  (B.  19)  as 

F    =  -  A-t-    ^B  -Y-;  n(y    .   )|  (B.  20) 

00  o  mm 

where    y    .      as  introduced  by  Griem,  Kolb  and  Shen  (Griem  et  al.  ,  1959) 
mm 

is  given  by 

4TTn  /       2\  ^  ^  I  2  \^ 

e  /  ehn   \  2  I         n  12, 

y    .    =    I    ,       I  =   -r  1-77  TtA   4C  .  B.21 

mm          3  m  V    kT   /  ^  \  3(^*1"^  4  )/ 

Consequently,  B  and  B      are  related  by  the  following  relation 


B    =  B  +Y+^ 


(B. 22) 


Comparing  Eqs.  (B.  19)  and  (B.  20)  one  notices  that  for  a  particular 

line  the  value  of  the  square  bracket  as  derived  here  depends  on  the 

quantum  number  n,  for  that  particular  state.     This  is  also  true  for  the 

k 

paper  of  Shen  and  Cooper,   1969,  who  consider  our  case  (b)  with  infinite 

limits  on  the  t' -integral.  Otherwise  the  constants  given  in  the  literature 

are  independent  of  n    because  the  lower  cutoff  parameter  is  usually 

2 

based  on  an  average  Stark  splitting.     If  we  set  nq-n'q'  =  n  /2,  which 
corresponds  approximately  to  the  average  Stark  splitting  and  also 
gives  the  results  for  the  Stark  shifted  component  of  Lyman-a,  we  have 


B    =  B-  0.  64        .  (B.  23) 

o 

2 

This  yields  directly  for  n^  -  n  /Z  the  B    values  corre'sponding  to  the 
B   values  in  Table  2. 
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The  following  constants  B    have  been  given  in  the  literature 

o 

Griem,  Kolb,  Shen,   1959;  (Eq.  29)   :  B    =  0 

o 

Griem,  Kolb,  Shen,   1962;  (Eq.  2)     :  B    =  1.0 

Griem,  1965 


Kepple,  Griem,  1968 


(neglecting  quadrupole  term)  :      B    =  0.  58 

o 


Shen,  Cooper,   1969;  B    =  0.  58 

o 

Recently  the  time  development  operator  (S-matrix)  has  been  evaluated 
for  Lyman-a  including  time  ordering  by  solving  the  differential  equations 
for  the  S-matrix  elements  (Bacon,  1969).,    Again  the  square  bracket 
depends  on  (nq-n'q')  and  the  average  value  B    =1,1  considering  only 
the  dipole  term.    It  should,  however,  be  stressed  that  one  should  not 
overinterpret  these  numbers  because  within  the  classical  path  approxi- 
mation there  is  always  some  uncertainty  about  the  "correct"  constant 
B   because  of  the  ambiguous  lower  cutoff.    This  is  due  to  the  fact 
that  the  classical  path  approximation  breaks  down  roughly  for  p<  ^ 
(for  details  see  Paper  I),    For  most  cases  this  has  no  significant  effect 
for  the  Stark  broadening  of  hydrogen  because  the  dynamic  broadening  is 
primarily  due  to  weak  collisions.    More  details  are  given  with  the 
discussion  at  the  end  of  Chap.  XII.    The  situation  is  quite  different  for 
the  broadening  of  ionized  lines  where  strong  collisions  are  very  im- 
portant and  where  the  uncertainty  of  the  classical  path  approximation 
accounts  for  part  of  the  still  existing  discrepancies  between  theory 
and  experiment,  which  are  large  compared  with  the  Stark  broadening 
of  hydrogen. 

So  far  we  have  considered  F  ,  which  is  the  basis  of  the  un- 

00 

modified  impact  theory.    In  order  to  extend  the  range  of  validity  beyond 
the  plasma  frequency  the  modified  impact  theory  introduces  the  Lewis 
cutoff  by  considering  only  those  collisions  for  which  the  duration  of  a 
collision,  which  is  typically  p/v,  is  smaller  than  the  time  of  interest 
being  typically  1/aud.    For  this  reason  the  modified  impact  theory 
introduces  an  upper  cutoff 
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Pmax=  ^^^<^'^a>^^  24) 


or 


It  should  be  noted  at  this  stage  that  in  the  following  relations  we  not 

only  have  to  replace  C  by  C/a  but  also  AU)_^  by  a   Auu_,  in  order  to  obtain 

R  R 

the  dependence  on  the  upper  cutoff  parameter  a.    Considering  the  us- 
ually applied  case  (b)  (Eq.  (B.  3))  we  have  to  evaluate  the  following 

integral  for  ^{u     >  I 
R 


duu^-y         dxx[cos(^-l]  ,B.25, 


u  AO)  T3  ^  / 

o     R  o/u 


where  the  lower  limit  on  the  u-integral  is  determined  by  the  condition 

u  /u<  1/auu„.    After  a  change  of  variables  and  a  partial  integration 
o  R 

one  obtains  similar  to  Eq.  (B.  9) 

^ZCAm^^ 

du  (B.26) 


1       r         -n'  (^^'">r\ 
:r  I        ue  cos  I  1  - 


2  Am 

R     u  auUt, 
o  R 


Expanding  the  cosine  again  and  performing  another  change  of  variables 
the  result  is 


2  2 

2       «  2k  °°      -u  Aw^z 

4 


2       «  2k  °°      -u  Aw^z 


1 

This  gives  us  then 


2  2k 


(B. 28) 
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Evaluating  the  exponential  integrals  E    for  small  arguments  only  we 
'nally  have 


L  2 


-^i(v^r)    +  ^1  ("o' 


+  I 


^2  3 

=  C    j-.n  AuOj^  +  i 


(B.29) 


''^■^'ch  gives  us  for  ^w  ^  >  1  the  log-dependence  of  the  $     -matrix  elements 

R  ab 

in  the  modified  impact  theory,    A  more  appropriate  way  for  applying  the 
Lewis  cutoff,  which  avoids  the  discontinuity  at  Aoj     =  1 ,   is  to  take  as 
an  upper  cutoff 


p  =  MIN  (D,  v/auu) 

max 


or  p        /D=MIN(1,  u/aoj^) 

max  R 


(B.  30) 


which  for  case  b   leads  to  the  following  integral. 


f  \  3  -uY  /  r  /2c\  1 

I    =1  du  u  e       f  dx  X     cos  I  1  -1 

Lj  J  I       \:.nj  J 

./  u^e- Y    ax  .  [cos  M  -:1 

A(«o  u  /u 

R  o 


(B.31) 
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These  integrals  are  identical  to 


duu  e      J  dxx 

V""  (B.32) 


0=  2  -V^^R 


-J     duu^-^  I       dxx    j^cos    ^-^^     -ij  . 


and  after  a  partial  integration  we  have 

00 

=  — —  I      du  u  e 


L  2 

AO) 


cos 


* '  o  R 


{B. 33) 


-u^r  /2c\  ,1 

+  i|     du  ue  cos  "1 


^'•"r  b 

Expanding  the  cosine  functions,  I  can  be  evaluated  in  a  similar 
manner  as  above  with  the  result 


_C_ 
2 


E^{A^j^)-2E^(u^Ac«j,) 


+  2C 


(2k)! 


k 


k=2 


2k-2 


(B.  34) 


E  (u  Auu„) 

2k-l      o  R 


00  2k-2 
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For  fSW-^  1    with  2C«1  and  b'<  1  the  latter  result  may  be  simplified 


to  give 


.2  r 


+  I 


(B. 35) 


b  b 

which  for  fsa^     ~*  0  reduces  down  to  I    =1     and  which  has  no  discontinuity 


at        =  1 

to 


Furthermore,  we  see  that  for  Aw 


R 


00  Eq.  (B.  34)  goes  over 
2 


u 


L,  Aw 


R 


2AaJ 


o  R 
e  cos 


R 


/2C\      e   /2C  \ 


(B. 36) 

which  does  not  lead  to  the  static  limit. 

Similar  results  can  be  obtained  for  case  a,  which  are  not 
included  because  they  are  no  longer  required.     The  derivation  for  case 
b  has  been  included,   in  order  to  obtain  consistent  relations  which  allow 
a  comparison  with  the  calculations  dons  in  this  paper.    The  results  for 
case    b   as  given  here  differ  slightly  from  the  results  in  the  literature 
which  also  vary  from  paper  to  paper  depending  on  the  average  matrix 
elements  used  and  on  what  lower  cutoff  and  average  velocity  has  been 
applied. 


102 


APPENDIX  C 

PROGRAM  FOR  CALCULATING  THE  FINAL  LINE  PROFILE  I(A(ju) 


In  the  following  a  complete  listing  of  the  program  is  given  which 
j  was  used  to  calculate  the  final  line  profile  I(Aa))  on  the  basis  of  the  unified 
!  theory  for  the  case  of  no  lower  state  interactions. 
1.     The  Fourier  transform  of  the  thermal  average 

The  complex  function  FOUTR  calculates  essentially  the  Fourier 
J  transform  of  the  thermal  average  as  defined  by 

FOUTR  =iAtUpi    (A(JUp,B,  n,  q      q  )  (C.l) 
K    u        K-  be 

i 

'  where 

DOM=  Auu„  =  (Auu-Aa).(  n,  q.  )  B)/uu  .  (C,2) 
R  X  b  p 

I  It  uses  the  Eqs.  (X.17)  and  (X.  22)  for  calculating  i  (A(i)„)  and  i  {f\(xi^) 
I  1       R  2  R 

respectively.     The  required  Bessel  functions  J    ,  J,  ,  Y    and  Y,  are 

o        1       o  1 

evaluated  by  the  subroutine  BSJYOl.     For  large  and  small  arguments 
i  ^^R   ^^^^^  relations  are  replaced  by  their  asymptotic  expansions  (X.  18), 
(X.  23)  and  (X.  26).     The  specifying  constants  p  ,  p  ,  b  ,  a    and  b 
(PI,  P2,   Bl,   A2,  B2)  are  set  in  the  function   AIIM  and   are  calculated 
once  for  all  the  Stark  components  in  the  main  program  STBRHY.  The 
function  FOUTR  can  be  replaced  by  another  short  function  FOUTR  listed 
at  the  end  of  the  program,  which  then   makes  the  program  calculate 
line  profiles  according  to  the  modified  impact  theory. 
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2.    Calculation  of  I(AijJ,  p) 

The  function  AIIM  calculates  I(AtJO,   B)  as  defined  by  Eq.  (XII.  6). 
It  establishes  first  of  all  the  matrix  of  the    jb-operator  according  to 
Eq.  (XII.  7)  and  calculates  the  array 


AMATR  (NBN„  NAN)  =  5         Auu„-(nq^mU(A(W     )Inqm)     (C.  3) 

q  q       R         b  op  a 

3,  D 


for  the  m  values  0  and  ±1.     The  required  3j- symbols  are  calculated 
once  in  the  main  program  and  their  values  are  stored  in  three  different 
arrays  according  to  the  following  definitions 

n-1       n-1  Y 


SAR  (NLA,  NQC)  =  (  Zl^  +  1  )^ 


and 


m 


I 


m  -q    m  +q 
c     c      c  c 


•m 


(C.4) 


i 


where    NLA  =  l    and  NQC  =  q 


SJQL  (MCT,  NBN,  NLA)  = 


/  n-J_         n-1  \ 


\ 


m-q  m+q 

a  a 


•m 


(C.5) 


/ 


SSJJ  (MCT,  NBN,  NAN) 

q_+q 

n+m-1-  — 

(-1) 


'b 


n-1 
2 

m-q 


n-1 
2 

m+q 


m 


/  n-1  n-_l 
2  2 


m-q^  m+q^ 


•m 


(  C.  6) 
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where  MCT  =  m  and   NAN  and  NBN   are  nurabering  indices  specifying 

the  q    and      .     The  matrix  AMATR  is  then  inverted  by  the  subroutine 
a  b  ' 

CGAUSSEL,  which  is  able  to  solve  systems  of  linear,  complex  equations 
by  Gaussian  elimination.    Multiplying  the  inverted  matrix  by  the  3j- 
symbols  according  to  Eq.  (XII.  6)  yields  finally  I{Aw,  6). 
3.     The  final  line  profile  I(Auu) 

The  main  program  STBRHY  calculates  the  final  line  profile  by 
performing  the  ion  field  average  according  to  Eq.  (II.  1).    It  first  of 
all  reads  in  the  ion  microfield  distribution  function  for  0  <  B  ^  30  in 
steps  of  0.  1,  which  has  been  calculated  in  a  separate  program  for  the 

particular  shielding  parameter  r  /D.    It  then  reads  in  the  density  n  , 

o  e 

the  temperature  T,  the  upper  principal  quantum  number  n,  the  wave- 
length X ,  the  average  value  of  the  static  ion  fields  8     ,  the  initial 
value  AUU)  the  logarithmic  stepwidth,  the  total  number  of  points,  a 
parameter  which  specifies  the  number  of  ion  field  integration  points 
and  finally  6  numbers,  which  specify  the  G  -function  and  hence  i  (aoj^) 
for  all  Stark  components  and  which  may  in  practically  all  cases  be  set 
to  zero.     These  6  numbers  are  obtained  from  the  thermal  average 
described  in  Appendix  A, 

The  program  then  calculates  the  constants  p^  ,  p^,  b^  ,   a^  and  b^ 
for  all  Stark  components  and  stores  them  in  the  array  FPAR.    p^  is 
calculated  according  to  Eq.  (X.  5).    p^   is  determined  on  the  basis  of 
Eq.  (B.  19)  where  the  constant  B  is  given  by 

B  =  0.27-  2K(-^)     .  (C.7) 
b 
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The  K-£unction  is  defined  in  Eq.  (B.  14)  and  is  calculated  for  a  lower 

cutoff  p     .    =  -X  +  -7  n  a     by  setting 
^  mm  O 


The  necessary  cosine  integral  is  calculated  by  the  function  COSINT. 


As  a  next  step  the  main  program  evaluates  all  the  required 
3j-symbols  by  means  of  the  function  S3  J,  which  in  turn  uses  the  function 
FCTRL,  to  calculate  all  the  necessary  factorials.    The  numbers  are 
stored  in  arrays  according  to  the  definitions  in  (C.  4)  ,  (C.  5)  and  (C.  6). 

In  performing  the  ion  field  average  the  microfield  distribution 

function  is  calculated  by  the  function  WFLD,  which  uses  a  5  point 

interpolation  on  the  values  read  in  initially  for  B<  30  and  otherwise 

uses  the  asymptotic  expansion  given  by  Hooper,  1968b.    As  a  function 

of  Auu  and  |3,  which  determine  the  shape  of  the  ion  field  integrand,  the 

ion  field  integral  is  subdivided  in  intervals,  which  are  integrated 

separately  by  means  of  Weddle's  rule  (subroutine  WEDDLE)  using  a 

convenient  change  of  variables  in  every  interval.     Furthermore,  the 

program  calculates  the  asymptotic  wing  expansion  according  to 

Eq.  (  XI.  12)  and  the  unified  theory  for   8=0  and  3  =  B^^  performing 

not  the  ion  field  integral.    All  three  values  are  normalized  with 

-5/2  . 

respect  to  the  asymptotic  -wing. 


106 


PROGRAM  STBRHY 

C 

C  PROGRAM  FOR  CALCULATING  THE  STARKBROADEN I NG  OF  HYDROGEN  ON  THE 

C  BASIS  OF   THE  UNIFIED  THEORY  FOR  NO  LOWER  STATE  INTERACTION 

DIMENSION   FF(1000)»   PFAC{6)»   SSJ(20),  STRONG(20) 

COMPLEX  FOUTR 

C0MM0N/FDAT/P1»P2»B1 »A2  »B2  »PPFF 
COMMON/ PSTR/NNN»NM1 »BET tFPAR ( 6 »20 ) 

COMMON/PSJD/SSJJ( 2»20»20) .   SJQL ( 2 »20 » 20 ) »  SAR(20»20) 
COMMON/ PFW/F I  ELD ( 301 ) 
FIELD(l)    =  0.0 

READ  100»    (FIELD(I)»    I   =  2,301) 
100     FORMAT  (6E12.4) 

120     READ  150,DEN»TEMP,NNN,ALAM»BAV»GIN»DGG»NTOT.NFAC» (PFAC{  I  ) »I  =  1  »6) 
150     FORMAT    (2E10.2,    15,   4F10,2,  2I5/6F10.5) 

IF    (EOF, 60)    577,  170 
170     PRINT    180,   DEN,   TEMP,   NNN»  ALAM 

180     FORMAT    (1H1»*     DENSITY  =  *E12,A*  TEMPERATURE   =  *E12.4, 

1  *       QUANTUMNUMBER   =*I2*       WAVELENGTH  =*F8.2*  ANGSTROM*// 

2  13X,»Pl*18X,*P2*18X,*Bl*18X,*A2*18X,*B2*17X,*STROWG*/ ) 
IF    (NNN.LE.20)    GO  TO  200 

PRINT  190 

190     FORMAT    (»  PROGRAM  NOT  EXECUTED  BECAUSE  N   IS  LARGER  THAN  20*) 

577     CALL  EXIT 

200     SDEN  =  SQRT{DEN) 

FAC  =   2064.936  *  TEMP  *  SQRT ( TEMP/DEN ) 

NMl   =   NNN  -  1 

NEVODD  =  M0D(NNN,2) 

AN  =  NNN 

AN2   =   AN  *  AN 

ANIM  =  0.5  *    (AN  -  1.) 

CFAC   =  4.5645E-7  *  AN  *  SDEN/TEMP 

DEBROG  =  2. 1027E-6/SQRT ( TEMP) 

RMIN  =  DEBROG  +  AN*AN*7 • 9376E-9 

DO  230  K  =   1,  NNN 

AK   =  K  —  1 

230     SSJ(K)=(aN2+( (-1. )**M0D(NMi+K,2) ) * ( AN2-2 . *AK*AK ) ) / ( 2 . *AN* ( AN2-1 • ) ) 
BET  =   5.6558E-5  *  AN  *  DEN**(l./6.) 

C 

C  ARRAY  FOR  G-FUNCTION  CONSTANTS 

ASY   =  0.0 
DO  270  <  =  1»NM1 
QC  =  K 

C  =  CFAC  *  QC 

PI   =  -1.671086  *  FAC  *  C  *  SQRT(C) 
BS     =   3.   *  AN  *  QC  *  DEBROG/RMIN 

STRONG (K)    =   0.269-2.* ( ( ( l.-COSF(BS) ) /BS+S INF ( BS ) ) /BS'COS iNT ( BS ) ) 

PPFF   =  -1.128379  ♦  FAC  *  C  *  C 

P2   =   PPFF  *    (STRONG(K)    -  2 . *LOGF ( 2 . *C ) ) 

FPAR(1,K)    =  PI 

FPAR{2»K)    =  P2 

FPAR{3»K)    =   0.5  *  (P2/P1)**2 
FIN   =  LOGF(AN*QC) 

FPAR(4,K)   =  P2  *   ( (PFAC(3)*FIN  +  PFAC(2))   *  FIN  +  PFAC(l)) 
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FPAR(5»K)    =    (PFAC(6)*FIN  +  PFAC(5))»FIN  +  PFAC(4) 

FPAR(6»K)    =  PPFF 
270     ASY   =   ASY  +  2.   *  FPAR(1,KJ    *  SSJ(K+l) 

PRINT   2A0 » ( { FPAR( K» I) =   1 ♦ 5 ) »STRONG ( I ) »    I   =  l.NMl) 
240     FORMAT  (6E20.4) 

PRINT   280»   FAC,   CFAC ,   BET*   ASY,   DEBROG»  NFAC 
?.8U     FORMAT(/*  FAC   =*E12.4,*  CFAC   =*E12.4,*  BET  =*E12.4,*  ASY  =*E12.4 

1  *  DEBR06  =*E12.4,*   INTEGRAT lONFACTOR  =*I2// 

2  5X*DOM*8X*DLAMi^8X*ITOT*6X*lHOLTS*8X*ASY*10X*WING*7X*WHOLTS*7X 

3  *WSTAT*7X*WWOO*8X*WWBB*8X*WWNG*/ ) 
AOLFAC   =  4.23538E-15  ^   SDEN  *  ALAM  *  ALAM 

C 

C  3JSYMB0L--ARRAYS  SS  J  J  (  MCT  ,  NBN  ,NAN  )    AND  S  JQL  (  McT  .NBN  ,NLa  ) 

DO  650  MCT  =1,2 
AMA  =  MCT  -  1 
NLIM   =  NNN  +   1   -  MCT 
NOB   =  -NLIM  -  1 
DO  450  NBN  =  1,NLIM 
NQB   =   NQB  +  2 
QB   =  NQB 

FMBl   =    (AMA  -  QB)    *  0.5 
FMB2   =    (AMA  +  QB)    *  0.5 
DO  325  NLA  =  1»NNN 
ALA  =   NLA  -  1 

325     SJQL(MCT,NBN,NLA)    =   S3 J ( AN IM ♦ AN IM » ALA » FMB 1 , FMB2 . -AMA ) 
450  CONTINUE 

DO  375  NBN  =  1,NLIM 

DO   350  NAN  =  1,NLIM 

AABB   =    (-1.  )^^*M0D(NAN+NBN.2  )    *  S  JQL  (  MCT  ,  NBN  ,  2  )    *  SJQL  (  MCT  .NaN  ,2  ) 
350     SSJJ(MCT,NBN,NAN)    =  AABB 
375  CONTINUE 
650  CONTINUE 

C 

C  3JSYMB0L-ARRAY  SAR(NLA»NQC) 

DO  780  NQC   =  l-NMl 
QC   =  NQC 

DO  720  NLA  =  1,NNN 
ALA  =   NLA  -  1 
FBB  =0, 

DO  680  NMC  =  1»NLA 

IF(NEV0DD.NE.M0D(NMC+NQC»2) )    GO  TO  680 
AMC   =   NMC  -  1 
FCF   =  2. 

IF{NMC.EQ.l)    FCF  =  1, 
FMCl   =   0.5  *    (AMC  -  QC) 
FMC2   =  0.5  *    (AMC  +  QC ) 

FBB  =   FBB  +  FCF  *    ( S3J ( ANIM, ANlM ♦ ALA tEMCl »FMC2 ,-AMC) ) **2 

680  CONTINUE 

720     SAR(NLA,NQC)    =   FBB  *   (2.*ALA  +  !•) 

780  CONTINUE 

C 

C  CALCULATION  OF   THE   lONFIELD  INTEGRAL 

NN6   =   6  *  NFAC 
ANN6   =  NN6 
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N12   =   12  *  NFAC 

AN12   =  N12 

N30   =   30  *  NFAC 

AN30   =  N30 

G  =  GIN  -  DGG 

DO  950  MM  =  1»NT0T 

G  =  G  +  DGG 

DOM  =   10.   **  G 

DLAM  =  ADLFAC  *  DOM 

WING52   =  -0.2992067103   *  AS Y/ ( SORT ( DOM )    *  DOM  *  DOM) 

FHOLTS  =  0. 

AWING  =  0. 

DO  815  NQC  =  1,NM1 

PI   =  FPAR(1»NQC) 

P2   =  FPAR(2»NQC) 

Bl   =  FPAR(3»NQC) 

A2   =  FPAR(4»NQC) 

B2   =  FPAR(5»NQC) 

PPFF   =  FPAR(6»NQC) 

AWING  =  AWING  +  SS J ( NQC  +  1 ) *A I  MAG ( FOUTR ( DOM )  ) *2 . / ( DOM*DOM ) 
QC  =  NQC 

BETFAC  =   BET   *  QC 
BCRIT   =  DOM/BETFAC 
815     FHOLTS   =   FHOLTS  +  SSJ ( NQC  +   1)    *  WFLD(BCRIT)    /  BETFAC 
AIRES   =  0. 

IF(DOM.GT. (-3.*P2 ) )    GO  TO  SAO 
BCRIT   =    (DOM  -  P2)/((AN  -   1 . ) *BET ) 
DB  =  BCRIT/ANN6 
B  =  0. 

DO  820  J  =  1»NN6 
B  =   B  +  DB 
820     FF(J)    =  AIIM(D0M,B)    *  WFLD(B) 

CALL  WEDDLE    ( DB ♦ NN6 ♦ FF . A  1 1  I » 0  .  ) 

AIRES  =  AI II 

DY   =  l./(BCRIT*AN30) 

Y  =  0. 

DO  830  J  =  1»N30 

Y  =   Y  +  DY 
B  =  l./Y 

830     FF(J)    =   B  »  B  *  AIIM(D0M»B)    *  WFLD(B) 

CALL  WEDDLE    ( DY ♦ N30 ♦ FF » A  I  I  I » 0  •  ) 

AIRES  =  AIRES  +  AIII 

GO  TO  980 
840     BCRCR  =  DOM/BET 

EPSPS  =  -P2/BET 

DO  957  NQ  =  1»NM1 

ANQ  =  NQ 

BCR   =  BCRCR/ANQ 

EPS   =  EPSPS/ANQ 

IF(NQ,EQ.l)    GO  TO  907 

SLl   =   l./(GAM  -  BCR) 

GO  TO  908 

907  SLl   =  0. 

908  SL2   =  l./EPS 
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SL3  =  l./(BCR  +  EPS) 
SL^  =   1,/(BCR  -  EPS) 

GAM  =  0.5  *    (SCR  -  EPS  +   (BCRCR  +  EPSPS)/(ANQ  +  !•)) 
IF(NQ.EQ.NM1 )   GAM  =  0.5  *   (BCR  -  EPS^ 
SL5  =   l./(BCR  -  GAM) 
CRIT  =  SL2  -  SL5 

Y  =  SLl 

IF(NQ.EQ.l)    GO  TO  913 
B  =  BCR  +  l./Y 

FA  =  AIIM(DQM,B)   *  WFLD(B)/(Y  *  Y) 
GO  TO  914 

913  FA  =  0. 

914  DY  =    (SL2  -  SL1)/ANN6 
DO  917  J  =  1»NN6 

Y  =  Y  +  DY 
Yl   =  l./Y 

B  =  BCR  +  Yl 
917     FF(J)    =  Yl  *  Yl  *  AIIM(DOM»B)   *  WFLD(B) 
CALL  WEDDLE   ( DY ^NNe »FF » A 1 1  I » FA ) 
AIRES   =  AIRES  +  AJII 

Y  =  SL3 
B  =  l./Y 

FA  =  B  *  B  ♦  AIIM(DOM»B)   *  WFLD(B) 
DY  =    (SL4  -  SL3)/ANN6 
DO  927   J   =  1»NN6 

Y  =  Y  +  DY 
B  =  l./Y 

927     FF(J)    =  B  *  B  »  AIIM(D0M,B)   *  WFLD(B) 
CALL  WEDDLE   ( DY .NN6 . FF » A  III » FA ) 
AIRES   =  AIRES  +  AI I  I 
IF(CRIT.LE.O. )   GO  TO  977 

Y  =  SL5 

B  =  BCR  -  l./Y 

FA   =  AIIM(D0M»8)    *  WFLD{B)/(Y  *  Y) 
DY   =  CRIT/AN12 
DO  937  J  =  1»N12 

Y  =   Y  +  DY 
Yl   =  l./Y 

B  =   BCR  -  Yl 
937     FF(J)    =  Yl  *  Yl  *  AIIM(DOM»B)   *  WFLD(B) 

CALL  WEDDLE    ( DY »N 12 .FF » A 1 1  I » FA ) 
957     AIRES  =  AIRES  +  AI  I  I 

IF(GAM.LT.   5.)    GO  TO  968 

Y  =  l./GAM 

FA   =  GAM  *  GAM  *  A I  I M ( DOM » GAM )    *  WFLD(GAM) 
DY  =    (0.2  -  Y)/AN12 
DO  967  J  =  1»N12 

Y  =   Y  +  DY 
B  =  l./Y 

967     FF(J)    =   B  *  B  *  AIIM(D0M»B)    *  WFLD(B^ 
CALL  WEDDLE    ( DY ♦N12 » FF » A  I  I  I , FA ) 
AIRES   =  AIRES  +  All  I 
SL4  =  0.2 
GO  TO  Q"77 
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966     SL4  =  1,/GAM 

977  B  =  0. 

[j  DB   =   l./(SL4  *  AN30) 

'  DO  947  J   =  1»N30 

B  =  B  +  DB 
947     FF(J)    =  AIIM{DOM»B)    *  WFLD(B) 

CALL  WEDDLE    ( DB  ♦NSO  ♦  FF » A I  I  I ♦ 0 . ) 
AIRES  =  AIRES  +  AI I  I 
980     WING   =  AIRES/WING52 

WISTAT   =  AlRES/FHOLTS 
WINHOL   =  FHOLTS/WING52 

WWOO   =    (AIIM(DOM»   0.)    +  FHOLTS ) / W I NG52 
WWBB   =    ( AIIM( DOM»BAV)    +  FHOL TS ) / W I NG52 
WWNG  =    (AWING  +  FHOLTS ) /WING52 
950     PRINT  978»D0M,DLAM,AIRES»FH0LTS»WING52»WING»WINH0L»WISTAT, 
2       WWOO*   WWBB»  WWNG 

978  FORMAT  (11E12.4) 
GO  TO  120 

END 

C 

FUNCTION  AIIM(DOM»B) 

I  C 

'C  CALCULATION  OF  I(DOM»B) 

COMPLEX  DFTR(2U) » AMATR ( 20 ♦ 40 ) »F0UTR,aREF 
COMMON/FDAT/Pl ,P2 ,31 »A2 ♦B2 ,PPFF 
C0MM0N/PSTR/NNN,NM1 ,BET»FPAR(6,20) 

C0MM0N/PSJD/S5JJ( 2»20,20) ♦   S JQL ( 2 ♦ 20 ♦ 20 ) ♦  SAR(20»20) 

AI  IM   =  0. 

DO  800  MCT   =  1,2 

AMCT   =  MCT 

NLIM   =  NNN  +   1   -  MCT 

NL22   =   2   *  NLIM 

NQB   =  -NLIM  -  1 

DO  750  NBN  =  1,NLIM 

NQB   =   NQB  +  2 

QB   =  NQB 

DOMRB  =  DOM  -  BET  »  QB  »  B 
DO  220  NQC  =  1,NM1 
PI   =  FPAR(1,NQC) 
P2   =  FPAR(2»NQC) 
Bl   =  FPAR(3,NQC) 
A2   =  FPAR(4,NQC) 
B2   =  FPAR(5»NQC) 
PPFF   =  FPAR(6,NQC) 
220     DFTR(NQC)    =  F0UTR(D0MRBJ 
DO  700  NAN   =  1»NLIM 
AREF   =  (0.,0,) 
DO  600  NQC   =  1,NM1 
FAA  =  0.0 

DO   500  NLA  =  1,NNN 
500     FAA  =   FAA  +  S JQL ( MCT ,NAN ,NLA )    *  S JQL ( MCT , NBN , NLA )    *  SAR(NLA,NQC) 
600     AREF   =   AREF  +  FAA  *  DFTR(NQC) 

AMATR{NBN,NAN  +  NLIM)    =  (0,,0.) 
700     AMATR(NBN,NAN)    =   6.2831853072  *    ( ( - 1 • ) **MOD ( NAN+NBN , 2 ) )    *  AREF 
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AMATR(NBN»NBN+NLIM)    =  (1,»0.) 
750     AMATR{NBN»NBN)    =  AMATR  (  NBN  ,NBN  )   -i-  DOMRB 
C  MATRIX  INVERSION 

CALL  CGAUSSEL( AMATR » 20 ♦NLIMjNLZZ^NRANK) 

DO  795  NBN  =  1»NLIM 

DO   793  NAN   =  1»NLIM 
793     AIIM   =  AIIM  +  AMCT*SSJJ ( MCT »NBN tNAN ) *A IMAG ( AMATR ( NBN »NAN+NL IM ) ) 
795  CONTINUE 
800  CONTINUE 

AIIM  =-AI IM  *  0.3183099 

RETURN 

END 

C 

C  THE   FOLLOWING  FUNCTION  FOUTR  MAY  BE  REPLACED  BY  THE   FUNCTION  FOUTf 

C  AT   THE  END  OF  THE  LISTING  TO  OBTAIN  THE  MODIFIED  IMPACTTHEORY 

FUNCTION  FOUTR  (DOM) 

C 

C  FOURIERTRANSFORM  OF  THERMAL  AVERAGE  FOR  UNIFIED  THEORY 

COMPLEX  FOUTR 

COMMON/FDAT/Pl ,P2 ,B1 ,A2»B2  »PPFF 
ARG  =  ABSF(DOM) 
Z   =   Bl   *  ARG 

IF    (Z.LE. 0.001)    GO  TO  600 
IF    (Z.LE. 40.)    GO   TO  300 

FACl   =  -0.2992067103   *  P 1 /( SORT ( ARG )    *  ARG  *  ARG) 
CC  =  FACl  *    ((1.   -  1.3125/Z)*0.625/Z  +  1.) 
SS   =   FACl   *    ((-1.-   1.3125/Z)*0.625/Z  +  1.) 

GO  TO  500 

300     CALL  BSJYOl    (Z»   AJO»   Y0»   AJl»  Yl) 
FACl   =   Y1/{2.*Z)    +  AJl   -  YO 
FAC2   =  AJO  +  Yl   -  AJ1/(2.*Z) 
CINE   =  COSF(Z) 
SINE   =  SINF(Z) 

CC  =  P2  *  Bl   *  Bl  *    (CINE  *  FACl   +  SINE  *  FAC2) 

SS  =  P2  *  Bl  *  Bl  *   (CINE  *  FAC2  -  SINE  *  FACl^ 

IF   (A2.EQ.0.)   GO  TO  500 

Z  =  B2  *  ARG 

IF    (Z.GT. 10. )    GO  TO  400 

CALL  BSJYOl    (Z»   AJO»   Y0»   AJl»  Yl) 

FACl  =  ( ( AJ1-Y0)*16.*Z-36.*AJ0-28.*Y1 )^Z+15.*YO-3.*AJl 
FAC2   =    ( ( AJO+Yl )*16.*Z-36.*Y0+28.*AJ1 )*Z-15.*AJ0-3.*Y1 

CINE  =  COSF(Z) 
SINE   =  SINF(Z) 

CC  =  CC  +  A2*B2*(CINE  *  FACl  +  SINE  *  FAC2)/6. 
SS  =  SS  +  A2*B2*(CINE  *  FAC2  -  SINE  *  FACl)/6. 
GO  TO  500 

400     FACl   =   0.1322319336  *  A2   *  B2  *  Z**{-3.5) 

CC   =  CC  +   FACl  *    (1.   -    (3.9375/Z  +  l.)*4.375/Z) 

SS   =  SS  -  FACl   *    (1.   -    (3.9375/Z  -  1,)*4.375/Z) 
500      IF    (DOM.LT.O.)    SS  =  -SS 

FOUTR  =  ARG  *  ARG  *  CMPLX ( -SS»CC ) 
RETURN 

600     SS   =    (P2*B1   -  A2)    *  DOM 

FOUTR   =  0.3183099  *  CMPLX ( SS »-P2 J 
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RETURN 
END 


C 
C 

c 


1 

2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 


SUBROUTINE  WEDDLE    ( DX ♦   N»   F»   A»   FO ) 


INTEGRATION  SUBROUTINE 
DIMENSION  F{N) 
A  =  0.0 
K  =   N  -  1 
DO   15    I   =    1»  6 
SUM  =  0#0 
DO  6  J  =    I »   K »  6 
SUM  =   SUM  +  F(J) 
GO  TO   (8»    lU,    12»    10.  8, 
5.0  * 


14)  ♦  I 


A  =  A  + 
GO  TO  15 
A  =  A 


SUM 


+  SUM 


200 


70 


74 
75 


GO  TO  15 

A  =  A  +  6.0  *  SUM 

GO  TO  15 

A  =  A  +  2.0  *  SUM 

CONTINUE 

A=0.3*DX*  <A+FO+F(N)) 

RETURN 

END 

FUNCTION  WFLD(B) 

CALCULATION  OF   THE   ION  MICROFIELD  DISTRIBUTION  FUNCTION  USING  A 
5P0INT    INTERPOLATION  FOR  THE   DATA  READ   INTO   THE  MAINPROGRAM 
COMMON/PFW/FIELD( 301 ) 
WFLD  =  0.0 

IF    (B.LE.30»0)    GO   TO  200 
SBS   =   l./(B  *  SQRT(B) ) 

WFLD  =    ((21.6  *  SBS  +  7.639)    *  SBS  +   1.496)   *  SBS/B 
RETURN 

IF  (B.LE.0.0)  RETURN 
J  =    (B  +  0.2 )    *  10.0 

L  =   J  -  1 

IF    (J.GT.2)    L  =  J  -  2 

IF    (J.GT.3)    L  =  J  -  3 

IF    (J.GT.300)    L   =  297 

LLL  =   L  +  4 

DO  75   K  =  L»LLL 

AK   =  K  -  1 

TERM  =  1.0 

DO  74  M  =  L»LLL 

IF    (K.EQ.M)    GO  TO  74 

AM  =  M  -  1 

TERM  =  TERM  *   (10. *B  -  AM)/(AK  -  AM) 
CONTINUE 

TERM  =   TERM  *  FIELD{K) 
WFLD  =  WFLD  +  TERM 
RETURN 
END 


c 

FUNCTION  S3J    (FJl,   FJ2»   FJ3»   FMl »   FM2 »   FM3  ) 

C 

C  CALCULATION  OF  3J-SYMB0L 

S3J=0#0 

IF(ABS(FM1  +  FM2  +  FM3 )    .GT.   0.001)    GO  TO  153 

FM3=FM1+FM2 

A=FJ2+FJ3+FMl+.O0  5 

B=FJl-FMl+.005 

C=-FjH-Fj2+FJ3+.005 

D=FJ3+FM3+«005 

E=FJl-FJ2-FM3+«005 

F=FJl-FJ2+FJ3+.005 

G=FJl+FJ2-FJ3+.005 

H=FJ1+FJ2+FJ3+1 .0+.005 

E2=FCTRL( B)*FCTRL(FJ1+FM1 )*FCTRL( FJ2-FM2 ) *FCTRL( FJ2+FM2) 
IF    (E2)  153»153»145 
145   E1=(FCTRL(C)*FCTRL(F)/FCTRL(H) ) *FCTRL ( G ) *FCTRL ( D ) *FCTRL ( F J3-FM3 ) 
IF    (El)    153»    153»  150 

150  E1=SQRT(E2)/SQRT(E1) 
I1=XMAX1F(0.0»  -E+0.01) 
I2=XMIN1F{A»  CD) 

IF    (12-11)    153,    151f  151 

151  DO   152  I=IliI2 
FI  =  I 

E2=FCTRL(FI J *FCTRL(C-FI )*FCTRL( D-FI ) / FCTRL ( A-F I ) 

152  S3J  =  S3J+(  n-l,0)**XMODF(  I»2)  )/E2)*FCTRL(B  +  FI )/FCTRL(E+FI ) 
U=ABS(FJ1+FM2+FM3 ) +0.001 

S3J=S3J*( (-1.0)**XMODF(XFIXF(U) »2) )/El 
FM3=-FM3 

153  RETURN 
END 

C 

FUNCTION  FCTRLCA) 

C 

C  CALCULATION  OF   FACTORIALS   REQUIRED  BY  FUNCTION  S3J 

DIMENSION   FCTI (20) 

C 

DATA    ( (FCTI  (  I  )  .  I  =  l»20)    = 1  .  0  ,  2  .  0  ♦  6 . 0  ♦  24 . 0  ♦  1 20 . 0  ♦  720 . 0 » 5040 . 0 ♦ 

1  40320.0  » 362 880.0 » 362 8 800.0 » 399 16800 •0*47900 1600.0 » 

2  622  7020  800 .0» 87 178291200.0 , 130  7674368000.0 , 

3  2.0922789888E13,    3  .  55687428096E 14 ,   6 . 402 373705 72 8E 1 5 ♦ 

4  1.2164510040883   E17»   2.432  9020081766  E18) 

C 

IF(A)  50»60»70 
50    IF( A.GE. (-0.1 )    )    GO  TO  60 

FCTRL   =  0.0 

RETURN 
60  FCTRL   =  1.0 

RETURN 
70    I   =   A  +  0.1 

IF    (    I    .EQ.   0)    GO  TO  60 

IF    (1-20)  140*140.130 
130  F=20.0 
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FCTRL=FCTI ( 20 ) 
DO   131  J=21»I 
F=F+1.0 

131  FCTRL=FCTRL*F 

GO  TO  150 
140   FCTRL  =  FCTI  {  I  ) 
150  RETURN 

END 

C 

SUBROUTINE  CGAUSSEL ( C » NRD ♦ NRR ,NCC » NSF ) 

C 

DIMENSION  C(NRD,NCC) »L( 128»2) 

TYPE  COMPLEX  C , DET » P ♦ D ♦ Q ♦ R 

DATA    (BITS  =    1755  4000   OOOO   0000  B) 

CALL  ROLLCALL (48HCGAUSSEL  6/5/68  l-BANK  BITS=2**-18  ) 
NR=NRR   $  NC=NCC 

IF(NC.LT,NR.OR.NR.GT,128.0R,NR.LE.O)   CALL  Q8QERR0R ( 0 ♦ 9HB AD  CALL.) 
C  INITIALIZE, 
NSF  =  0 

NRM=NR-1   S  NRP=NR+1   S  D=(1.»0.)    $  LSD=1 

DO  1  KR=1»NR 

L(KR.l)=KR 

1  L(KR.2)=0 

CALL  Q9EXUN(EXUN) 

IF(NR.EQ.l)    GO  TO  42 
C  ELIMINATION  PHASE, 

DO  41  KP=1»NRM 

KPP=KP+1   S   PM=0.    $  MPN=0 
C  SEARCH  COLUMN  KP   FROM  DIAGONAL  DOWN,   FOR  MAX  PIVOT. 

DO  2  KR=KP»NR 

LKR=L (KR» 1 ) 

PT=CABS(C(LKR»KP) ) 

IF(PT.LE.PM>    GO  TO  2 

PM=PT   S  MPN=KR  $  LMP=LKR 

2  CONTINUE 

C  IF  MAX  PIVOT   IS  ZERO,  MATRIX   IS  SINGULAR. 

IFCMPN.EQ.O)    GO  TO  9 
NSF=NSF+1 

IF(MPN.EQ.KP)    GO  TO  3 
C  NEW  ROW  NUMBER  KP  HAS  MAX  PIVOT. 

LSD=-LSD 

L(KP»2 ) =L(MPN,1 )=L{KP,1 ) 
L(KP,l)=LMP 

C  ROW  OPERATIONS  TO  ZERO  COLUMN  KP  BELOW  DIAGONAL. 

3  MKP=L(KP,1) 
P=C(MKP,KP)    S  D=D*P 
DO  41  KR=KPP,NR 
MKR=L(KR) 

Q=C {MKR,KP ) /P 

IF(REAL(Q) .EQ.O..AND.AIMAG(Q) .EQ.O. )    GO  TO  41 

C  SUBTRACT  Q  *  PIVOT  ROW  FROM  ROW  KR. 

DO  4  LC=KPP»NC 
R=Q*C{MKP,LC) 
C(MKR,LC> =C(MKR»LC)-R 
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4  IF(CABS(C(MKR»LC) ) .LT,CABS(R)*BITS)    C ( MKR » LC » = ( 0 . , 0 . ) 

41  CONTINUE 

C  LOWER  RIGHT  HAND  CORNER. 

42  LNR=L(NR,1)    S  P=C{LNR»NR) 

IF(REAL(P) .EQ.O..AND.AIMAG(P) .EQ.O, )    GO  TO  9 

NSF=NSF+1 

D=D*P*LSD 

IF{NR,EQ.NC)    GO  TO  8 
C  BACK  SOLUTION  PHASE. 

DO  61  MC=NRP»NC 
C(LNR.MC)=C(LNR,MC) /P 
IF(NR.EQ.l)    GO   TO  61 
DO  6  LL=1»NRM 

KR=NR-LL  S  MR=L{KR,1)    $  KRP=KR+1 

DO   5  MS=KRP»NR 

LMS=L{MS»1) 

R=C (MR»MS )*C ( LMSjMC) 

C(MR.MC)=C(MR,MC)-R 

5  IF(CABS(C(MR»MC) ) .LT.CABS(K)*BITS)    C ( MR » MC ) = ( 0 . , 0 . ) 

6  C(MR»MC)=C(MR,MC) /C(MR»KR) 
61  CONTINUE 

C  SHUFFLE  SOLUTION  ROWS  BACK   TO  NATURAL  ORDER. 

DO  71  LL=1»NRM 
KR=NR-LL 
MKR=L(KR»2) 
IF(MKR.EQ.O'    GO  TO  71 
MKP=L(KR,1) 
DO  7  LC=NRP»NC 
Q=C(MKR»LC) 
C(MKR,LC)=C(MKP,LC) 

7  C{MKP,LC)=Q 
71  CONTINUE 

C  NORMAL  AND  SINGULAR  RETURNS.      GOOD  SOLUTION  COULD  HAVE  D=0. 

8  C( 1 ,1 )=D  $  GO  TO  91 

9  C(1.1)=(0.»0.) 

91  CALL  S9FAULT{ EXUN)    $  RETURN 

END 

C 

SUBROUTINE  BSJYOl    (X»   AJO»    YO ♦   AJ1»  YD 

C 

C  CALCULATION  OF   THE  BESSEL  FUNCTIONS  JO ♦   YO,   Jl,   AND  Yl   FOR  AN 

C  ARGUMENT  X 

DIMENSION  A(7)»   B{7)»   C(7),   D(7),   E(7),   F(7).   G(7)»  H(7) 

DATA    ((All),    I    =   1,7)    =     0.00021,   -0.0039444,  0.0444479, 
1       -0.3163866,    1.2656208,   -2.2499997,  1.0) 

DATA    ((B(I)»    I    =   1,7)    =  -0.00024846,   0.00427916,  -0.04261214, 
1        0.25300117,   -0.74350384,   0.60559366,  0.36746691) 

DATA    (IC(I)»    I    =   1,7)    =     0.00014476,   -0.00072805,  0'.00137237, 
1       -0.00009512,   -0.0055274,   -0.00000077,  0.79788456) 

DATA    ((D(I)»    I    =    1,7)    =     0.00013558,   -0.00029333,  -0.00054125, 
1        0.00262573,   -0.00003954,   -0.04166397,  -0,78539816) 

DATA    ((E(I)»    I    =   1,7)    =     0.00001109,   -0.00031761,  0.00443319, 
1       -0.03954289,   0.21093573,   -0.56249985,  0.5) 

DATA    ((F(I)»    I    =    1,7)    =     0.0027873,   -0.0400976,  0.3123951, 
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1       -1,3164827»   2.16827U9*   0.2212091,  -0.6366198) 

DATA    ((G(I)»    I   =   1»7)    =  -0, 00020033*   0,00113653»  -0.00249511» 
1       0.0O017105»   0.01659667*   0,00000156»  0.79788456) 

DATA    ((H(I)»    I   =   1»7)    =  -0.00029166,   0.00079824,  0.00074348, 
1       -0.00637879,   0.0000565,   0.12499612,  -2.35619449) 

AX   =  ABSF(X) 

IF    (AX.GT.Q.O)    GO   TO  10 

AJO  =1.0 

YO  =  -l.E+030 

AJl   =  0.0 

Yl   =  -l,E+030 

RETURN 

10        IF    IAX.GT.3.0)    GO  TO  50 

XX   =    (AX/3.0)    **  2 

AJO  =  All) 

YO  =  Bd) 

AJl   =  Ed) 

Yl   =   F(  1  ) 

DO  20  M  =  2»7 

AJO  =  AJO  *  XX  +  A(M) 

YO   =  YO  *  XX  +  B(M) 

AJl   =  AJl  *  XX  +  E(M) 
20       Yl   =   Yl  *   XX  +  F(M) 

AJl   =   AJl  *  X 

ALF   =   0.6366197724  *  LOGF(0.5  *  AX ) 

YO  =  YO  +  ALF  *  AJO 

Yl   =  Yl/X  +  ALF  *  AJl 

RETURN 
50       X3   =  3.0/AX 

FO   =  C(  1) 

THO   =   D( 1  ) 

Fl   =  G( 1) 

THl   =  H(l ) 

DO  60  M  =  2»7 

FO   =   FO  *   X3   +  C(M) 

THO  =   THO  *  X3  +  D(M) 

Fl   =   Fl   ♦   X3   +  G(M) 
60       THl   =   THl   *  X3  +  H(M) 

THO   =   THO  +  AX 

THl   =   THl   +  AX 

XS   =  l./SQRT{AX) 

AJO=   XS  *  FO  *  COSF(THO) 

YO   =   XS  *   FO  *  SINF(THO) 

AJ1=   XS  ♦   Fl   *  COSF(THl) 

Yl   =  XS  *   Fl  *  SINF(THI) 

RETURN 

END 

C 

FUNCTION  COSINT(X) 

C 

C  CALCULATION  OF   THE  COSINE  INTEGRAL 

TYPE  DOUBLE  Y2 , PROD , SUM , PT , DK 
IF(X.LE.O.)   GO  TO  50 
X2   =  X  *  X 
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IF(X.GT.20.^    GO  TO  30 

Y2   =  DBLE(X2) 
PROD  =  -Y2  *  0.5 
SUM  =  PROD  *  0.5 
DO  10  K  =  2»50 
DK  =  2  *  K 

PROD  =  -PROD  *  Y2/(DK*(DK  -  1.)) 

SUM  =  SUM  +  PROD/DK 

PT   =  ABSIPROD  *  1.D+10> 

IF(ABS(SUM) .GT.PT)   GO  TO  20 
10  CONTINUE 
20   SS   =  SNGL(SUM) 

COSINT   =  SS     +  0.5772156649  +  LOGF(X> 

RETURN 
30  FA  =  1. 

FB  =  1. 

PO  =  1. 

X2  =  1./X2 

DO  40   K  =  1»10 

AK   =   2   *  K 

PO  =  -PO  *  AK  ♦  X2 

FA  =   FA  +  PO 

PO  =  PO  ^^    (AK  +  1.) 

FB  =  FB  +  PO 

PA  =  ABS(PO  *  I.E+IO) 

IF(PA.LE.FB)    GO   TO  45 
40  CONTINUE 
45   FX   =  FA/X 

GX   =   FB  *  X2 

COSINT   =   FX  ♦  SIN(X)    -  GX  *  COS(X) 
RETURN 
50  WRITE    (61»100)  X 
100  FORMAT    (*  X  LESS  OR  EQUAL  TO  ZERO,   X  =  *E17.9) 

RETURN 
END 

FUNCTION  FOUTR  (DOM) 

FOURIERTRANSFORM  OF  THERMAL  AVERAGE  FOR  MODIFIED  IMPACTTHEORY 
COMPLEX  FOUTR 

COMMON/FDAT/Pl ,P2,B1 ,A2 ,B2 ♦PPFF 

ARG  =  ABSF(DOM) 

CC  =  P2  ^  ^„ 

IF    (ARG. GT. 1.22474)    CC   =  P2  -  2 . *LOGF ( ARG/ 1 . 22474 ) *PPFF 
FOUTR   =  0.3183099  *  CMPLX ( 0 . » -CC ) 

•  RETURN 
END 
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Figure  1.      Schematic  picture  of  the  collision  sphere  showing  the 
Debye  sphere,  a  strong  collision  sphere  and  a  straight  line  classical 
path  trajectory. 
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Figure  2.    The  thermal  average  F  of  the  time  development  operator 
normalized  with  respect  to  the  static,   small  interaction  time  asymp- 
tote F    as  a  function  of  the  normalized  time  s  =*(ju    •  t.    The  two  curves 
are  obtained  with  two  different  lower  ^utoff  paran?eters  in  the  p- 

integral,  n     .    =  0     and  n     .    =  ^  +  n  a  . 

mm  mm  o 
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F(S) 
Fo  (S) 


Log,o  S 

Figure  3.     The  thermal  average  F  of  the  time  development  operator 
normalized  with  respect  to  the  static,   small  interaction  time  asymp- 
tote F    as  a  function  of  the  normalized  time  s  =  %  t.     The  two  sets  of 
o  p 

curves  are  obtained  with  two  different  lower  cutoff  parameters  in  the 

2 

p -integral,  p     .    =  0  and  p     .    - +  n  a    .     The  three  different  curves 
mm  mm  o 

in  every  set  correspond  to  different  Stark  components  characterized  by 

the  quantum  number  n^  =  nq-n'q'. 
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Log,o  Aoip  — - 

Figure  4.    The  Fourier  transform  of  the  thermal  average  i,  normalized 
with  respect  to  the  static,  large  frequency  limit  i    as  a  function  of  the 
normalized  frequency  Aujp  =  (Atu-Auu.«  B)/uj'  . 
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n,^  =  90  16  3 


Figure  5.    The  Fourier  transform  of  the  thermal  average  i,  normalized 
with  respect  to  the  static,  large  frequency  limit  i    as  a  function  of  the 
normalized  frequency  Aui)-^  =  (Aw-AU).*  P •    The   three  different  curves 
correspond  to  different      Stark  corn.pone:&s  characterized  by  the  quan- 
tum number  n,  =  nq~n'q'. 
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Figure  13.  The  final  line  profiles  for  the  density  and  temperature 
parameters  of  Vidal,  1965. 
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Figure  14.    Plot  of  the  electron  density  values  as  a  function  of  the 
principal  quantum  number  n,  which  have  been  evaluated  by  Vidal, 
1965,  under  the  assumption  that  the  AX~^/^  -wing  revealed  by  the 
experiment  is  identical  with  the  asymptotic  Holtsmark  /^,X"^/^-  wing 
for  electrons  and  ions. 
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Figure  15.  The  Balmer  line  profiles  for  the  density  and  temperature 
parameters  of  Vidal,  1965,  as  calculated  by  the  indicated  authors. 
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Figure  16.    The  Lyman-p  profile  for  two  different  strong  collisions 

cutoffs  p     .  . 

mm 
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THE  NATIONAL  ECONOMIC  GOAL 

Sustained  maximum  growth  in  a  free 
market  economy,  without  inflation, 
under  conditions  of  full  employment 
and  equal  opportunity 


THE  DEPARTMENT  OF  COMMERCE 

The  historic  mission  of  the  Department 
is  "to  foster,  promote  and  develop  the 
foreign  and  domestic  commerce"  of  the 
United  States.  This  has  evolved,  as  a 
result  of  legislative  and  administrative 
additions,  to  encompass  broadly  the  re- 
sponsibility to  foster,  serve  and  promote 
the  nation's  economic  development  and 
technological  advancement.  The  Depart- 
ment seeks  to  fulfill  this  mission  through 
these  activities: 


MISSION  AND 
FUNCTIONS 
OF  THE 

DEPARTMENT  OF 
COMMERCE 


"to  foster,  serve  and 
promote  the  nation's 
economic  development 
and  technological 
advancement" 


Participating  with 
other  government 
agencies  in  the 
creation  of  national 
policy,  through  the 
President's  Cabinet 
and  its  subdivisions. 

•  Cabinet  Committee 
on  Economic  Policy 

•  Urban  Affairs 
Council 

•Environmental 
Quality  Council 


Promoting  progressive 
business  policies  and 
growth. 

•  Business  and 
Defense  Services 
Administration 

•  Office  of  Field 
Services 


Assisting  states, 
communities  and 
individuals  toward 
economic  progress. 

•  Economic 
Development 
Administration 

•  Regional  Planning 
Commissions 

•  Office  of  Minority 
Business  Enterprise 


NOTE:  This  schematic  Is  neither  an  organization  chart  nor  a 
progrann  outline  for  budget  purposes.  It  is  a  general  statement 
of  the  Department's  mission  in  relation  to  the  national  goal 
of  economic  development. 
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Strengthening 
the  international 
economic  position 
of  the  United 
States. 

•  Bureau  of 
International 
Commerce 

•  Office  of  Foreign 
Commercial 
Services 

•  Office  of  Foreign 
Direct  Investments 

•  United  States 
Travel  Service 

•  Maritime 
Administration 


Assuring  effective 
use  and  growth  of  the 
nation's  scientific 
and  technical 
resources. 

•  Environmental 
Science  Services 
Administration 

•  Patent  Office 

•  National  Bureau  of 
Standards 

•  Office  of 
Telecommunications 

•  Office  of  State 
Technical  Services 


Acquiring,  analyzing 
and  disseminating 
information  concern- 
ing the  nation  and 
the  economy  to  help 
achieve  increased 
social  and  economic 
benefit. 

•  Bureau  of 
the  Census 

•  Office  of  Business 
Economics 
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NBS  TECHNICAL  PUBLICATIONS 


PERIODICALS 

JOURNAL  OF  RESEARCH  reports  National 
Bureau  of  Standards  research  and  development  in 
physics,  mathematics,  chemistry,  and  engineering. 
Comprehensive  scientific  papers  give  complete  details 
of  the  work,  including  laboratory  data,  experimental 
procedures,  and  theoretical  and  mathematical  analy- 
ses. Illustrated  with  photographs,  drawings,  and 
charts. 

Published  in  three  sections,  available  separately: 

•  Physics  and  Chemistry 

Papers  of  interest  primarily  to  scientists  working  in 
these  fields.  This  section  covers  a  broad  range  of 
physical  and  chemical  research,  with  major  emphasis 
on  standards  of  physical  measurement,  fundamental 
constants,  and  properties  of  matter.  Issued  six  times 
a  year.  Annual  subscription:  Domestic,  $9.50;  for- 
eigii,  $11.75*. 

•  Mathematical  Sciences 

Studies  and  compilations  designed  mainly  for  the 
mathematician  and  theoretical  physicist.  Topics  in 
mathematical  statistics,  theory  of  experiment  design, 
numerical  analysis,  theoretical  physics  and  chemis- 
try, logical  design  and  programming  of  computers 
and  computer  systems.  Short  numerical  tables. 
Issued  quarterly.  Annual  subscription:  Domestic, 
$5.00;  foreign,  $6.25*. 

•  Engineering  and  Instrumentation 

Reporting  results  of  interest  chiefly  to  the  engineer 
and  the  applied  scientist.  This  section  includes  many 
of  the  new  developments  in  instrumentation  resulting 
from  the  Bureau's  work  in  physical  measurement, 
data  processing,  and  development  of  test  methods. 
It  will  also  cover  some  of  the  work  in  acoustics, 
applied  mechanics,  building  research,  and  cryogenic 
engineering.  Issued  quarterly.  Annual  subscription: 
Domestic,  $5.00;  foreign,  $6.25*. 

TECHNICAL  NEWS  BULLETIN 

The  best  single  source  of  information  concerning  the 
Bureau's  research,  developmental,  cooperative  and 
publication  activities,  this  monthly  publication  is 
designed  for  the  industry-oriented  individual  whose 
daily  work  involves  intimate  contact  with  science  and 
technology — for  engineers,  chemists,  physicists,  re- 
search managers,  product-development  managers,  and 
company  executives.  Annual  subscription:  Domestic, 
$3.00;  foreign,  $4.00*. 

•  Difference  in  price  is  due  to  extra  cost  of  foreign  mailinff. 

Order  NBS  publications  from: 


NONPERIODICALS 

Applied  Mathematics  Series.  Mathematical  tables, 
manuals,  and  studies. 

Building  Science  Series.  Research  results,  test 
methods,  and  performance  criteria  of  building  ma- 
terials, components,  systems,  and  structures. 

Handbooks.  Recommended  codes  of  engineering 
and  industrial  practice  (including  safety  codes)  de- 
veloped in  cooperation  with  interested  industries, 
professional  organizations,  and  regulatory  bodies. 

Special  Publications.  Proceedings  of  NBS  confer- 
ences, bibliographies,  annual  reports,  wall  charts, 
pamphlets,  etc. 

Monographs.  Major  contributions  to  the  technical 
literature  on  various  subjects  related  to  the  Bureau's 
scientific  and  technical  activities. 

National  Standard   Reference   Data  Series. 

NSRDS  provides  quantitive  data  on  the  physical 
and  chemical  properties  of  materials,  compiled  from 
the  world's  literature  and  critically  evaluated. 

Product  Standards.  Provide  requirements  for  sizes, 
types,  quality  and  methods  for  testing  various  indus- 
trial products.  These  standards  are  developed  coopera- 
tively with  interested  Government  and  industry  groups 
and  provide  the  basis  for  common  understanding  of 
product  characteristics  for  both  buyers  and  sellers. 
Their  use  is  voluntary. 

Technical  Notes.  This  series  consists  of  communi- 
cations and  reports  (covering  both  other  agency  and 
NBS-sponsored  work)  of  limited  or  transitory  interest. 

Federal  Information  Processing  Standards  Pub- 
lications. This  series  is  the  official  publication  within 
the  Federal  Government  for  information  on  standards 
adopted  and  promulgated  under  the  Public  Law 
89-306,  and  Bureau  of  the  Budget  Circular  A-86 
entitled.  Standardization  of  Data  Elements  and  Codes 
in  Data  Systems. 

CLEARINGHOUSE 

The  Clearinghouse  for  Federal  Scientific  and 
Technical  Information,  operated  by  NBS,  supplies 
unclassified  information  related  to  Government-gen- 
erated science  and  technology  in  defense,  space, 
atomic  energy,  and  other  national  programs.  For 
further  information  on  Clearinghouse  services,  write: 

Clearinghouse 

U.S.  Department  of  Commerce 
Springfield,  Virginia  22151 
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